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Abstract 


The  effect  of  Mach  number  on  the  plane  mixing  layer  has  been  investigated 
by  means  of  linear  stability  theory  and  two-  and  three-dimensional  direct  numerical 
simulations  of  the  compressible  Navier-Stokes  equations.  The  objective  was  to  iden¬ 
tify  the  effects  of  compressibility  on  a  building-block  fluid  flow,  with  applications 
to  supersonic  mixing  and  combustion. 

Results  from  linear  stability  theory  show  that  the  amplification  rate  is  reduced 
as  Mach  number  is  increased.  Above  a  convective  Mach  number  of  0.6  it  is  found 
that  three-dimensional  waves  are  more  amplified  than  two-dimensional  waves  and 
a  simple  relation  is  found  to  give  the  orientation  of  the  most  amplified  waves.  It  is 
also  shown  that  the  linear  stability  theory  can  be  used  to  predict  the  mixing  layer 
growth  rate  as  a  function  of  velocity  ratio,  density  ratio  and  Mach  number. 

Two-dimensional  simulations  show  a  strong  reduction  in  growth  rate  of  the  two- 
dimensional  motion  as  Mach  number  is  increased,  with  more  elongated  structures 
forming  at  high  Mach  numbers.  Shock  waves  are  observed  in  two-dimensional  sim¬ 
ulations  above  a  convective  Mach  number  of  0.7.  The  supersonic  modes  of  insta¬ 
bility,  which  are  the  only  two-dimensional  unstable  modes  at  high  Mach  numbers, 
are  shown  to  be  radiating  and  vortical,  but  have  very  low  growth  rates. 

Three-dimensional  simulations  with  random  initial  conditions  confirm  the  linear 
stability  result  that  oblique  waves  become  the  most  amplified  waves  at  high  Mach 
numbers,  with  no  evidence  for  any  other  modes  of  instability.  Simulations  beginning 
with  a  two-dimensional  wave  and  a  pair  of  equal  and  opposite  oblique  waves  show 
a  change  in  the  evolved  large-scale  structure  as  Mach  number  is  increased.  Above 
a  convective  Mach  number  of  0.6  the  oblique  modes  have  most  of  the  energy  in  the 
developed  structure,  and  above  a  convective  Mach  number  of  1  the  two-dimensional 
instability  wave  has  little  effect  on  flow  structure.  Similar  organized  structure  was 
found  in  a  simulation  with  random  initial  conditions.  No  shock  waves  were  found 
in  the  three-dimensional  simulations,  even  at  convective  Mach  numbers  above  1. 
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CHAPTER  1 


Introduction 


1.1  Motivation 

Free  shear  layers  are  of  fundamental  importance  in  the  study  of  turbulence.  They 
are  found  in  many  situations  in  nature  ( e.g .  atmospheric  flows,  volcanic  eruptions, 
stellar  jets)  and  in  industrial  applications  (e.g.  gas  turbine  combustor,  airfoil  wake, 
rocket  exhaust).  Detailed  understanding  of  the  physics  of  free  shear  layers  is  essen¬ 
tial  for  the  development  of  new  turbulence  and  mixing  models.  Improved  models 
of  mixing  in  free  shear  layers  will  lead  to  a  better  capability  for  the  prediction  of 
chemical  reactions  and  control  of  pollutant  emissions,  for  example  from  oil  and  gas 
burners  in  power  generation  plant. 

Progress  in  space  research  is  dependent  on  developing  more  efficient  propulsion 
systems,  and  vehicles  capable  of  carrying  a  higher  payload  into  orbit.  Future  fully 
re-usable  space  vehicles  have  been  proposed,  which  would  be  capable  of  taking  off 
from  conventional  runways  and  attaining  earth  orbit.  Such  vehicles  would  use  air 
breathing  engines,  and  in  the  range  of  flight  Mach  numbers  from  5  to  20  it  appears 
(Swithenbank  et  at.  [1989])  that  the  best  efficiency  is  obtained  in  a  supersonic 
combustion  ram  jet  (scramjet).  In  such  an  engine  the  heat  addition  takes  place 
at  supersonic  speeds  and  the  air  velocity  throughout  the  engine  is  approximately 
equal  to  the  flight  velocity.  The  limiting  process  in  these  engines  is  the  time  taken 
to  mix  the  fuel  and  oxidizer,  which  must  occur  within  the  combustion  chamber  for 
the  heat  release  to  be  of  value  in  generating  thrust.  Mixing  occurs  in  supersonic 
free  shear  layers  within  the  combustor.  A  general  understanding  of  the  physics  of 
compressible  mixing  may  suggest  methods  of  reducing  the  mixing  time  and  making 
scramjets  more  efficient. 

The  plane  mixing  layer  is  a  simple  prototype  of  a  free  shear  layer,  consisting  of  two 
streams  of  fluid  with  unequal  velocities,  and  in  the  compressible  case  often  unequal 
densities.  The  low  speed  version  has  been  extensively  studied  in  the  laboratory, 
and  compressible  experiments  are  now  being  performed.  The  flow  is  amenable  to 
solution  by  direct  numerical  simulation,  and  has  been  selected  as  the  basic  flow  for 
this  study. 
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1.2  Survey  of  Previous  Work 


Much  previous  research  has  been  directed  towards  the  incompressible  mixing 
layer,  but  the  effect  of  compressibility  has  received  relatively  little  attention.  For 
completeness,  the  literature  for  both  high  and  low  speed  mixing  layers  is  reviewed 
here. 

1.2.1  Experiments 

Early  mixing  layer  research  mapped  out  the  time-averaged  character  of  the  flow 
(Liepmann  and  Laufer  [1947]),  and  identified  the  strong  dependence  of  the  down¬ 
stream  development  of  the  flow  on  upstream  effects  (Bradshaw  [1966]).  A  ma¬ 
jor  change  in  thinking  about  the  mixing  layer  occurred  when  Brown  and  Roshko 
[1974]  observed  that  large-scale  nearly  two-dimensional  structures,  which  had  pre¬ 
viously  been  associated  with  a  transition  phenomena,  persisted  in  the  flow  at  high 
Reynolds  numbers,  when  the  mixing  layer  was  statistically  self-similar.  The  large 
scale  structure  has  been  found  by  many  succeeding  researchers,  including  Dimotakis 
and  Brown  [1976]  at  even  higher  Reynolds  numbers.  Oster  and  Wygnansky  [1982] 
found  that  two-dimensional  disturbances  applied  at  the  splitter  plate  were  able  to 
control  the  appearance  of  the  large-scale  structures  downstream  in  the  mixing  layer. 
They  identified  regions  where  growth  was  enhanced  or  retarded  by  the  effects  of  the 
forcing. 

Three-dimensional  structure  was  observed  in  the  original  work  of  Brown  and 
Roshko  [1974],  showing  up  as  streamwise  streaks  in  the  braid  region  between  suc¬ 
cessive  rollers.  This  secondary  structure  was  investigated  in  detail  by  Bernal  and 
Roshko  [1986],  who  showed  that  it  consisted  of  counter-rotating  streamwise  vor¬ 
tices  in  the  braids,  the  ends  of  which  became  wrapped  around  the  neighboring 
large  rollers.  The  counter-rotating  vortices  move  fluid  between  them  up  or  down 
which  shows  up  as  a  mushroom  shaped  structure  in  the  scalar  field,  shown  clearly 
in  Bernal’s  pictures.  The  streamwise  vortices  were  initially  fixed  in  space,  devel¬ 
oping  from  small  rig-dependent  disturbances  in  the  upstream  flow  field.  Later  in 
the  mixing  layer  development  they  appear  to  move  around  in  space,  since  long 
time-exposure  pictures  did  not  show  their  presence. 

Identification  of  dominant  structures  in  the  incompressible  mixing  layer  has  led 
to  the  development  of  new  models  capable  of  predicting  the  behavior  of  the  flow. 
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Dimotakis  [1986]  proposed  a  model  for  entrainment,  based  on  the  geometry  of  the 
primary  two-dimensional  motion,  which  successfully  predicts  the  experimentally- 
observed  phenomenon  of  asymmetric  entrainment  of  fluid,  with  more  fluid  from 
the  high-speed  side  than  the  low-speed  side  being  entrained  into  the  mixing  layer. 
An  important  new  model  for  mixing  and  chemical  reactions  has  been  developed  by 
Broadwell  and  Breidenthal  [1982],  and  for  reactions  with  finite-rate  chemistry  by 
Broadwell  and  Mungal  [1988].  The  model  uses  observations  from  experiment  to 
identify  various  fluid  states  in  the  mixing  layer:  (i)  unmixed  fluid,  (ii)  fluid  in  the 
structures,  mixed  at  the  entrainment  ratio,  and  (iii)  fluid  in  strained  laminar  diffu¬ 
sion  layers  between  the  two  free  streams  ( e.g .  in  the  braid  region).  Mixing  in  each 
region  has  its  own  characteristic  behavior  as  a  function  of  Reynolds  and  Schmidt 
numbers.  The  Broadwell-Breidenthal-Mungal  model  correctly  predicts  many  trends 
in  experiments  with  chemical  reactions  (Mungal  and  Dimotakis  [1984],  Breidenthal 
[1981],  Koochesfahani  and  Dimotakis  [1986],  Mungal  and  Frieler  [1988]),  and  has 
been  used  to  predict  the  effects  of  forcing  on  mixing  (Sandham  et  al.  [1988]). 

The  effect  of  compressibility  on  the  plane  mixing  layer  was  first  investigated 
in  the  1960’s,  for  the  mixing  layer  between  one  high-speed  stream,  and  another 
stream  at  rest.  The  data,  compiled  by  Birch  &  Eggers  [1973],  showed  a  reduction 
in  the  growth  rate  of  the  mixing  layer  as  the  Mach  number  was  increased.  Brown  & 
Roshko  [1974]  found  that  the  density  ratio  alone  could  not  account  for  the  reduction 
in  growth  rate,  implying  that  a  true  compressibility  effect  was  being  observed. 

Renewed  interest  in  compressible  mixing  in  the  1980’s  led  to  experiments  by 
Bogdanoff  [1983],  and  by  Papamoschou  <fc  Roshko  [1986,  1988],  in  which  high¬ 
speed  mixing  layers  with  various  velocity  and  density  ratios  were  investigated.  Both 
sets  of  researchers  proposed  a  parameter,  called  the  convective  Mach  number  by 
Papamoschou  &  Roshko,  which  seemed  to  collapse  all  the  available  growth  rate 
data  onto  one  curve,  showing  reduced  growth  rate  as  the  convective  Mach  number 
was  increased.  The  reasoning  behind  the  convective  Mach  number  concept  can  be 
found  in  Papamoschou  and  Roshko  [1988].  It  is  based  on  the  existance  of  organized 
large-scale  structure  in  the  compressible  mixing  layer.  If  one  looks  at  the  flow  in  a 
reference  frame  convecting  with  the  large  structures,  then  the  Mach  number  of  the 
free-streams  is  an  intrinsic  Mach  number  for  the  flow.  Respectively  for  each  stream, 


we  have: 


Mei  = 


(1.1) 
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where  U*  is  the  convective  velocity  of  the  structures,  U*  and  t/|  are  the  free-stream 
velocities,  and  cj  and  c£  the  sound  speeds.  The  superscript  *  denotes  a  dimensional 
quantity.  In  the  incompressible  mixing  layer  in  a  reference  frame  moving  at  Uc 
there  is  a  stable  stagnation  point  in  the  braid  region.  If  the  existence  of  a  similar 
point  is  assumed  for  the  compressible  layer,  and  that  this  point  is  reached  by  an 
isentropic  process  from  the  free-streams,  then  an  expression  for  U*  can  be  derived 
(Papamoschou  and  Roshko  [1988]).  For  gases  with  71  =  72  it  is  found  that  Mc\  = 
Mc 2  and: 


u;=m±m 


'’i+'i 


(1.2) 


We  can  also  eliminate  U*  from  equation  (l.l),  and  define  the  Mach  number  Mc  as: 


~  U2 
1 +  c2 


(1.3) 


More  recently,  Papamoschou  [1989]  has  attempted  to  measure  the  convective  ve¬ 
locities  of  the  large  scale  structures  directly  from  experimental  Schlieren  images. 
He  finds  disagreement  between  experiment  and  theory,  and  has  proposed  an  alter¬ 
native  two-dimensional  large-scale  structure,  in  which  shock  waves  are  allowed  on 
one  side  of  the  mixing  layer,  breaking  the  assumption  of  isentropic  flow  in  the  above 
derivation. 

Recent  experiments  by  Samimy  et  al.  [1989]  show  a  reduction  in  turbulence  levels 
as  Mach  number  is  increased.  Ongoing  flow-visualization  work  at  Stanford  (Clemens 
et  al.  [1989])  shows  little  evidence  for  organized  two-dimensional  motion  at  a  con¬ 
vective  Mach  number  of  0.6. 

1.2.2  Linear  Stability  Theory 

Numerical  solutions  of  the  linearized  equations  first  began  to  appear  in  the  early 
1960’s.  Applications  to  the  incompressible  mixing  layer  were  presented  by  Michalke 
[1965a, b,c],  and  to  the  compressible  mixing  layer  by  Lessen  et  al.  [1965,1966],  and 
by  Gropengiesser  [1970].  Earlier  analytical  work  had  revealed  the  instability  of  a 
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velocity  profile  with  an  inflection  point  (Rayleigh  [1880])  and  the  preference  for 
amplification  of  two-dimensional  waves  in  incompressible  flows  (Squire  [1933]). 

The  valuable  contributions  of  Gropengiesser  [1970]  were  largely  overlooked  at 
the  time.  He  found  the  second  mode  of  instability,  previously  observed  by  Lessen 
et  al.  [1966],  which  keeps  the  mixing  layer  unstable  in  two  dimensions  at  high 
Mach  number.  This  instability  was  subsequently  rediscovered  by  Blumen,  Drazin 
et  al.  (1970,  1975,  1977].  Gropengiesser  used  spatial  stability  theory  and  used 
a  solution  of  the  compressible  laminar  boundary-layer  equations  as  the  base  flow, 
instead  of  a  simple  hyperbolic  tangent  profile.  He  noted  the  high  amplification  rate 
of  three-dimensional  waves  at  high  Mach  number,  as  found  for  the  compressible 
wall  boundary-layer  by  Mack  (see  e.g.  Mack  [1984]). 

When  there  are  walls  present  in  the  flow,  or  in  wakes  with  an  embedded  sub¬ 
sonic  region  relative  to  the  free-stream,  there  can  be  additional  modes  of  instability 
present.  These  ‘acoustic’  modes  were  found  by  Mack  [1989]  in  wall  boundary  layers 
and  in  near  wakes.  Additional  modes  were  found  for  the  confined  mixing  layer  by 
Greenough  et  al.  [1989],  referred  to  by  them  as  ‘wall  modes’,  which  may  be  the 
same  kind  of  instability. 

The  important  effects  of  the  mean  velocity  profile  were  investigated  by  Monkewitz 
and  Huerre  [1982].  They  found  that  only  the  amplification  rate  computed  by  spatial 
theory  for  the  Blasius  mixing  layer  velocity  profile  showed  growth  rate  proportional 
to  A  =  (U£  —  U*)/{U*  +  U^),  as  found  in  experiments.  Morkovin  [1988]  makes  the 
point  that  only  the  results  from  spatial  stability  analysis  based  on  a  mean  profile 
satisfying  the  boundary-layer  equations  can  be  compared  with  experiments.  Earlier 
attempts  to  transform  results  from  the  temporal  theory  gave  poor  agreement  with 
the  measurements. 

Huerre  and  Monkewitz  [1985]  introduced  the  notion  of  convective  and  absolute 
instabilities  for  free  shear  layers,  which  affects  the  choice  of  temporal  versus  spatial 
theory.  If  the  flow  is  convectively  unstable  the  linear  theory  for  spatially  growing 
disturbances  is  applicable,  whereas  the  temporal  theory  should  be  used  for  abso¬ 
lutely  unstable  flows.  They  found  that  the  transition  from  absolute  to  convective 
instability  for  the  mixing  layer,  unfortunately  based  on  a  hyperbolic  tangent  veloc¬ 
ity  profile,  occurred  at  X  =  1.315,  t.e.  for  flows  with  a  significant  backflow  on  one 
side. 
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Recent  work  by  Ragab  and  Wu  [1988]  shows  that,  as  with  incompressible  free 
shear  layers,  the  compressible  mixing  layer  has  a  basically  inviscid,  inflectional 
instability  and  the  effect  of  viscosity  is  only  to  damp  the  disturbances.  They  also 
found  that  non-parallel  effects  are  negligible  in  compressible  mixing  layers. 

1.2.3  Secondary  Stability  Theory 

Secondary  stability  theory  was  developed  to  allow  a  further  step  into  transition, 
beyond  the  primary  instability,  to  be  studied.  The  theory  assumes  that  the  primary 
instability  has  developed,  modifying  the  basic  flow  field.  A  new  eigenvalue  stability 
problem  is  set  up,  in  which  the  base  flow  and  eigenfunction  are  dependent  upon  both 
the  streamwise  and  the  normal  location.  This  approach  has  proven  very  successful 
for  wall  boundary-layer  and  channel  flows,  predicting  the  appearance  of  both  K 
(Klebanoff)  and  H  (Herbert)  type  breakdown  towards  turbulence  (Herbert  [1983], 
Herbert  and  Bodonyi  [1989]),  which  have  been  observed  in  experiment  and  in  direct 
numerical  simulation  (Singer  et  al.  [1987]). 

The  major  work  in  this  field  for  the  mixing  layer  was  performed  by  Pierrehumbert 
and  Widnall  [1982].  They  assumed  a  base  flow  consisting  of  the  hyperbolic  tangent 
mixing  layer  mean  flow,  with  superposed  Stuart  [1967]  vortices,  which  are  steady 
solutions  to  the  incompressible  Navier-Stokes  equations.  Two  classes  of  instability 
modes  were  found:  fundamental  and  subharmonic.  The  fundamental  modes  have 
the  same  wavelength  in  the  streamwise  direction  as  the  vortex  spacing,  and  the  sub¬ 
harmonic  modes  have  twice  the  wavelength.  Two  fundamental  modes  were  found, 
corresponding  to  vortex  core  deformations  -  a  core  ‘bulging’  mode  (spanwise  sym¬ 
metric),  and  a  core  ‘translative’  mode  (spanwise  antisymmetric).  The  translative 
mode  was  the  more  unstable,  and  the  wavelength  of  the  most  (rapidly)  amplified 
instability  was  roughly  equal  to  the  spanwise  spacing  of  streamwise  vortices  found 
in  the  experiments  of  Bernal  and  Roshko  [1986].  The  most  amplified  subharmonic 
wave  was  the  two-dimensional  subharmonic,  corresponding  to  the  2D  pairing  pro¬ 
cess  most  commonly  observed  in  experiments.  Pierrehumbert  and  Widnall  did  find 
another,  three-dimensional,  subharmonic  mode.  This  produced  a  helical  pairing, 
and  may  possibly  have  been  observed  in  the  flow  visualizations  of  Chandrsuda  et 
al.  [1978]. 

The  only  work  in  secondary  stability  for  the  compressible  mixing  layer  was  car¬ 
ried  out  by  Ragab  and  Wu  [1989].  They  used  a  base  mixing  layer  profile,  and 
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superimposed  the  neutral  mode  of  the  two-dimensional  instability.  They  studied 
the  subharmonic  instability,  and  found  that,  above  a  convective  Mach  number  of 
0.4,  the  helical  pairing  mode  was  more  amplified  than  the  2D  pairing  mode.  This 
is  a  further  indication  that  the  dominant  instability  modes  at  high  Mach  number 
are  three-dimensional  modes. 

1.2.4  Numerical  Simulations 

Techniques  for  time-accurate  numerical  simulation  of  free  shear  layers  can  be 
divided  into  three  categories  (i)  vortex  dynamics  calculations,  (ii)  large-eddy  simu¬ 
lations  and  (iii)  direct  numerical  simulations.  The  first  ( e.g .  Ashurst  and  Meiberg 
[1988])  assumes  inviscid,  incompressible  flow.  The  Biot-Savart  rule  for  vortex  induc¬ 
tion  assumes  an  instantaneous  transfer  of  information  in  the  flowfield,  which  does 
not  happen  in  compressible  flow,  where  the  speed  of  sound  is  finite.  The  method 
of  large-eddy  simulation  (LES)  requires  a  model  for  the  smallest  scales  of  turbu¬ 
lence.  The  method  of  direct  numerical  simulation  (DNS)  can  produce  spatially  and 
temporally  accurate  solutions  of  the  full  Navier-Stokes  equations  with  no  modeling, 
when  care  is  taken  to  choose  the  flow  parameters  (t.g.  Reynolds  number,  Schmidt 
number)  in  order  to  fully  resolve  all  the  scales  of  motion.  Usually  the  requirement 
for  spatial  resolution  of  a  wide  range  of  scales  necessitates  the  use  of  spectral  or 
very  high  order  finite-difference  numerical  methods. 

Two  types  of  mixing  layer  problem  can  be  tackled  numerically.  The  spatially- 
developing  mixing  layer  computations  (figure  1.1)  use  the  same  reference  frame  as 
the  experiments,  with  inflow  at  one  end  of  the  computational  box  and  outflow  at  the 
other  end.  These  inflow/outflow  boundaries  require  special  treatment  -  given  the 
convective  nature  of  the  instability,  the  computed  mixing  layer  will  be  dominated 
by  the  upstream  forcing,  as  specified  by  the  inflow  boundary  condition.  The  outflow 
boundary  must  allow  all  structures  to  smoothly  leave  the  computational  domain, 
without  reflection  of  waves  back  into  the  simulation.  The  alternative  computations 
are  for  the  time-developing  mixing  layer  (figure  1.2).  Here  the  computational  do¬ 
main  is  fixed  in  a  reference  frame  moving  with  the  structures,  and  the  flow  then 
develops  in  time,  rather  than  in  space.  Periodic  boundary  conditions  are  enforced 
in  the  streamwise  direction  and  the  flow  develops  from  a  specified  initial  condition. 
Time-developing  simulations  permit  a  more  efficient  use  of  computational  resources, 
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and  the  periodic  boundary  conditions  in  the  streamwise  and  spanwise  direction  are 
tailor-made  for  highly  accurate  Fourier  methods. 

The  earliest  mixing  layer  computations  were  performed  with  LES  for  the  time- 
developing  mixing  layer  (Mansour  et  al.  [I978j,  Cain  et  al.  [1981]).  The  first  di¬ 
rect  numerical  simulations  were  presented  by  Riley  and  Metcalfe  [1980]  for  a  time- 
developing  layer  developing  from  random  initial  conditions,  using  methods  devel¬ 
oped  by  Orszag  and  Pao  [1974].  Recent  work  (Metcalfe  et  al.  [1987],  Rogers  and 
Moser  [1989])  has  shown  that  the  experimentally  observed  phenomena  of  primary 
roll-up  and  secondary  streamwise  vortices  leading  to  mushroom  shaped  structures 
in  the  braids  can  be  fully  realized  in  numerical  simulations. 

The  spatially-developing  mixing  layer  was  simulated  by  Lowery  and  Reynolds 
[1986].  Lowery  used  inflow  boundary  conditions  consisting  of  the  mean  flow,  with 
eigenfunctions  of  the  fundamental  and  two  subharmonic  frequencies  (from  linear 
stability  analysis)  superimposed.  This  resulted  in  the  generation  of  a  forced  mix¬ 
ing  layer  and  he  was  able  to  show  that  the  asymmetry  of  entrainment,  observed 
in  experiments,  could  be  captured  in  a  spatially-developing  computation.  Lowery 
also  performed  three-dimensional  computations  in  which  streamwise  vortices  were 
present  in  the  inflow  field.  These  developed  the  characteristic  mushroom  structure 
in  the  braid  region.  Follow-up  work  by  Sandham  and  Reynolds  [1989],  using  Low¬ 
ery’s  2D  code,  showed  that  the  large  asymmetry  of  entrainment,  observed  in  the 
initial  development  of  the  mixing  layer  (Koochesfahani  and  Dimotakis  [1986])  could 
be  traced  to  the  effect  of  the  wake  of  the  splitter  plate,  upstream  of  the  development 
of  the  mixing  layer.  It  was  also  shown  that  a  random-walk,  applied  to  the  phase  of 
the  forcing  eigenfunctions,  could  be  used  at  the  inflow  to  simulate  a  more  natural 
mixing  layer,  with  linear  growth  rates  and  more  nearly  self-similar  statistics. 

The  approach  of  Corcos  et  al.  (Corcos  and  Sherman  [1984],  Corcos  and  Lin  [1984] 
and  Lin  and  Corcos  [1984]),  has  been  to  directly  simulate  simple  deterministic  model 
flows  of  the  mixing  layer,  and  extract  physical  information  from  these.  Part  I  of 
their  study  deals  with  the  two-dimensional  roll-up,  Part  II  with  the  development 
of  the  secondary  instability,  and  Part  III  with  the  effect  of  the  staining  field  from 
the  primary  roll-up  on  the  development  of  streamwise  vorticity.  They  found  that 
the  plane  strain  field,  produced  in  the  braid  region  by  the  roll-up  of  the  primary 
instability,  acted  to  ‘collapse’  the  streamwise  vorticity  associated  with  the  secondary 
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instability  of  the  flow,  into  circular  streamwise  vortices,  which  then  generated  the 
mushroom  structures  in  the  scalar  field  in  the  braid. 

Direct  numerical  simulations  of  the  compressible  mixing  layer  have  only  recently 
appeared.  Both  the  two-dimensional  simulations  in  the  current  work,  and  those 
of  Lele  [1988,  1989]  show  the  reduction  in  growth  rate  of  the  mixing  layer  as  the 
convective  Mach  number  is  increased,  and  the  appearance  of  weak  embedded  shock 
waves  in  the  flow  for  convective  Mach  numbers  above  0.7.  Lele  also  shows  that 
simulations  started  from  low  level  random  perturbations  evolve  first  into  the  usual 
primary  roll-up,  with  wavelength  getting  longer  as  the  Mach  number  is  increased, 
as  predicted  by  the  linear  stability  analysis.  Other  work  is  limited,  but  Soestrisno 
et  al.  [1988]  presented  two-dimensional  simulations  of  the  time-developing  mixing 
layer,  and  Eberhardt  et  al.  [1988]  have  presented  simulations  of  the  ‘wall  mode’ 
of  instability  for  the  confined  mixing  layer.  The  simulations  show  that  this  mode 
tends  to  kink  the  mixing  interface,  but  does  not  lead  to  a  roll-up,  and  probably 
would  not  contribute  to  enhanced  mixing  at  high  Mach  number. 


1.3  Objectives  and  Overview 

The  objective  of  this  work  is  to  secure  a  fundamental  understanding  of  the  effect 
of  compressibility  on  the  development  of  the  plane  mixing  layer.  Several  important 
questions  emerge  from  the  previous  work  in  this  area.  (1)  Why  does  the  mixing 
layer  grow  more  slowly  at  higher  Mach  numbers?  (2)  Is  the  convective  Mach  number 
a  good  parameter  to  describe  the  complete  behavior  of  the  flow?  (3)  What  is  the 
structure  of  the  largest  scales  in  the  mixing  layer  at  high  Mach  number,  and  what 
implications,  if  any,  does  this  have  for  mixing  at  supersonic  speeds? 

The  approach  of  this  work  is  numerical.  First,  in  Chapter  2,  the  linear  stability 
problem  is  formulated  and  solved  for  a  wide  variety  of  mixing  layers.  Spatial  stabil¬ 
ity  analysis  is  used  to  compare  with  experiments  and  temporal  stability  to  obtain 
eigenfunctions  for  use  in  the  direct  numerical  simulations  in  later  chapters.  Inter¬ 
esting  features  emerge  in  their  own  right  from  the  stability  computations.  First, 
it  is  shown  that  the  amplification  rate  of  the  most  unstable  mode  from  the  linear 
stability  analysis  can  be  used  to  correctly  predict  the  growth  rate  of  the  developed 
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mixing  layer.  This  suggests  that  linear  processes  may  be  very  important,  even  in 
the  developed  mixing  layer. 

Secondly,  it  is  found  that  three-dimensional  waves  are  far  more  amplified  at  high 
Mach  numbers  than  two-dimensional  waves,  implying  that  the  developed  structure 
at  high  Mach  number  is  highly  three-dimensional. 

Results  from  linear  theory  are  used  to  provide  eigenfunctions  for  initial  conditions 
for  the  numerical  simulations.  Linear  theory  is  also  used  to  indicate  the  key  Mach 
numbers  at  which  the  instability  characteristics  of  the  flow  change,  where  changes  in 
physics  may  be  observed,  and  to  suggest  Reynolds  numbers  to  use  in  the  simulations, 
low  enough  to  resolve  the  flow,  but  high  enough  to  capture  the  inviscid  nature  of 
the  instabilities. 

The  numerical  methods,  used  for  the  direct  numerical  simulations,  are  presented 
in  Chapter  3.  The  three-dimensional  code  is  spectral  (Fourier)  in  the  periodic 
directions,  and  high-order  compact  finite  difference  in  the  normal  direction. 

In  Chapter  4  two-dimensional  simulations  are  presented,  which  illustrate  the  re¬ 
duction  in  growth  rate  as  the  Mach  number  is  increased.  Consideration  of  the 
compressible  vorticity  equation  shows  how  dilatational  and  baroclinic  effects  can 
explain  the  stabilization  of  the  two-dimensional  instability,  as  compressibility  be¬ 
comes  marked.  At  higher  Mach  numbers  the  two-dimensional  simulations  show  the 
appearance  of  weak  shock  waves,  embedded  around  the  vortices. 

The  key  three-dimensional  effects  are  presented  in  Chapter  5.  At  low  Mach 
number  the  modes  found  by  Pierrehumbert  and  Widnall  [1982]  are  simulated,  by 
carefully  choosing  the  phasing  of  a  pair  of  oblique  instability  waves  relative  to  the 
fundamental  2D  wave.  The  effect  of  increasing  Mach  number  on  the  structure  of  the 
mixing  layer  is  investigated  by  running  three  simulations  at  Mach  numbers  where 
(i)  the  2D  mode  is  dominant,  (ii)  both  2D  and  3D  waves  are  approximately  equally 
amplified,  and  (iii)  high  Mach  number  where  the  3D  waves  are  dominant.  The  latter 
two  simulations  identify  new  structures  in  the  mixing  layer,  which  are  highly  three- 
dimensional.  A  model  structure  for  high  Mach  number  flow  (Afc  >  1)  is  developed, 
based  on  the  non-linear  development  of  two  equal  and  opposite  oblique  waves. 
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The  main  contributions  of  this  work  are: 

•  mixing  layer  growth  rate  can  be  predicted  by  linear  stability  analysis.  The 
mixing  layer  growth  rate  is  found  to  be  directly  proportional  to  the  linear  am¬ 
plification  rate  of  the  most  amplified  spatial  instability  wave,  using  a  solution 
to  the  boundary-layer  equations  for  the  base  velocity  and  temperature  profiles. 

•  the  most  amplified  instability  wave  becomes  an  oblique  wave  above  a  convec¬ 
tive  Mach  number  of  0.6.  At  higher  Mach  numbers  the  most  amplified  wave 
becomes  more  oblique,  and  the  relation  Mc  cos  6  =  0.6  was  found  to  predict 
the  orientation  of  the  most  amplified  waves  at  high  convective  Mach  numbers. 

•  the  non-linear  growth  rate  of  the  two-dimensional  instability  wave,  which  dom¬ 
inates  the  incompressible  mixing  layer,  is  reduced  as  Mach  number  is  increased. 

•  if  the  flow  is  forced  to  be  two-dimensional  then  shock  waves  develop  for  con¬ 
vective  Mach  numbers  above  0.7. 

•  three-dimensional  instability  waves  at  high  Mach  number  grow  strongly  in  the 
non-linear  region  of  roll-up,  as  well  as  in  the  linear  regime.  The  developed 
structure  is  found  to  change  as  Mach  number  is  increased,  with  less  spanwise 
coherence  and  strong  streamwise  vorticity  at  higher  convective  Mach  numbers. 

•  the  structure  that  develops  at  high  Mach  number  from  a  pair  of  equal  and 
opposite  oblique  instability  waves  consists  of  a  pair  of  hairpin-like  vortical 
structures  which  are  split  in  a  peak-valley  manner  in  the  streamwise  direction. 
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CHAPTER  2 


Linear  Stability  Theory 


In  this  chapter  the  linearized  theory  describing  the  growth  of  small  disturbances 
in  the  compressible  mixing  layer  is  considered.  The  linearized  equations  are  solved 
to  find  the  most  amplified  wave  for  given  profiles  of  velocity  and  density.  It  will 
be  shown  how  linear  stability  theory  can  be  used  to  predict  the  growth  rates  of 
mixing  layers,  giving  good  agreement  with  experiments  for  the  effects  of  velocity 
ratio,  density  ratio  and  Mach  number.  Three-dimensional  waves  are  found  to  be 
important  in  the  compressible  mixing  layer  at  high  Mach  number.  The  structure  of 
the  eigensolutions  gives  important  information  about  the  structure  and  growth  of 
the  vortices  which  develop  out  of  the  linear  instability. 

In  the  following  sections  the  inviscid  equations  are  used  and  parallel  flow  is  as¬ 
sumed.  Spatial  theory  is  used  when  comparisons  with  experiments  are  desired. 
Temporal  theory  is  used  when  eigenfunctions  are  desired  as  inputs  to  direct  numer¬ 
ical  simulations  of  the  time-developing  mixing  layer. 

2.1  Numerical  Solution  Schemes 


The  compressible  boundary-layer  equations  are  solved  using  a  shooting  technique 
to  obtain  the  mean  flow.  The  linear  disturbance  equations  are  then  solved  using 
a  shooting  procedure.  The  methods  are  basically  from  Gropengiesser  [1970],  but 
have  been  extended  to  allow  computation  of  both  temporal  and  spatial  instability 
characteristics  of  a  variety  of  planar  free  shear  layers. 

2.1.1  Solution  for  the  Mean  Flow 


The  boundary-layer  equations  for  steady  two-dimensional  flow  of  a  compressible 
perfect  gas  with  zero  streamwise  pressure  gradient  are  (White  [1974]): 


W  djpV) 

dx*  dy* 


(2.1) 


.  ,9u‘  ,  ,du‘  9  /  ,du‘\ 

’“a?  +  ',,’d?-=a?V'‘  aP) 
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*  %dh*  ,  tdh*  a  / 

j  u - 1-  p  v  - = - ( 

dx*  p  dy *  dy *  \ 


dy* 


d 

(H*  dh*s 

dy* 

\Pr  dy*) 

(2.3) 


The  streamwise  direction  is  x*  and  the  normal  direction  y* .  Velocity  components 
in  these  directions  are  u*  and  v*  respectively.  The  density  is  denoted  by  p* ,  the 
viscosity  by  /z*  and  the  enthalpy  per  unit  mass  by  h* .  The  superscript  *  represents 
a  dimensional  quantity  and  the  Prandtl  number  is  defined  by  Pr  =  c*p.* /k* ,  where 
k*  is  the  conductivity. 


The  perfect  gas  law  is 


p*  =  p'R'T' 


(2.4) 


where  R *  is  the  gas  constant  and  p*  is  the  pressure.  Constant  specific  heats  are 
assumed,  and  we  write  dh*  =  CpdT* ,  where  is  the  specific  heat  at  constant 
pressure. 

Non-dimensionalization  of  the  above  equations  is  obtained  by  dividing  the  di¬ 
mensional  quantity  by  the  corresponding  dimensional  quantity  on  the  high-speed 
(y  >  0)  side  of  the  mixing  layer  (subscript  1).  The  new  dimensionless  variable  has 
no  superscript.  For  example: 


u  = 


(2.5) 


The  reference  length-scale  is  for  the  time  being  an  arbitrary  constant  /*,  and  the 
reference  time-scale  is  l* /U^.  Non-dimensional  parameters  are  the  Reynolds  number 
Rei  =  the  Mach  number  of  the  high-speed  side  of  the  mixing  layer 

Mi  =  U^/c^,  where  the  sound  speed  is  denoted  by  c*,  and  the  ratio  of  specific 
heats  7  =  c*/c*. 

The  boundary-layer  equations  in  dimensionless  form  become: 


d{pu)  d{pv) 
dx  dy 


=  0 


(2.6) 


du  du  1  d  (  du\ 

l,urI+l‘v-^  =  Rrl^{llai) 


(2.7) 
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and 


pu 


dT 

dx 


dT 


1  d  dT  Mfo  - 1)  /du\2 
ReiPrdy^dy  Rc\  dy / 


(2.8) 


pT  =  1  (2.9) 

since  pressure  has  been  assumed  uniform. 

The  first  step  is  to  derive  a  relation  between  temperature  and  velocity.  The  pro¬ 
cedure,  described  in  White  [I974j,  is  to  search  for  a  solution  T  =  T(u).  Substituting 
into  the  energy  equation  (2.8),  and  using  for  example 

dT  _  dTdu 
dy  du  dy 


we  have 


dT  i  du  du 

*r(p“^  +  ""ai 


i  a 

ReiPr  dy 


(2.11) 


If  the  Prandtl  number  is  assumed  unity  then  the  left  hand  side  must  be  zero  by  the 
momentum  equation  (2.7),  and  we  have  the  equation 


Mfb-U  (2.12) 

This  can  be  integrated  twice,  subject  to  the  boundary  conditions  that  T  =  1  when 
u  =  1  and  r  =  r2  when  u  =  f/2 ,  giving: 

r = M*2 -r1  (“i1 + -  “2  -  +  +  (f^tl  (2-13) 

This  is  the  general  T  —  u  relation,  referred  to  as  a  Crocco-Busemann  relation  after 
the  original  developers. 

The  next  step  is  to  find  a  solution  to  the  continuity  and  momentum  equations. 
A  stream  function  rj),  which  automatically  satisfies  continuity,  is  defined  by: 


(2.14) 
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The  momentum  equation  now  becomes: 


dip  du  dip  du  1  d  du 
dy  dx  dx  dy  Rei  dy*1  dy 


(2.15) 


We  now  introduce  Howarth’s  transform  (also  known  as  the  Illingworth-Stewartson 
transform,  see  Schlichting  [1979]): 


Using  this,  together  with  equation  (2.9),  we  can  rewrite  (2.15)  as: 

dip  d^ip  dip  d^ip  1  d  n  d2ip 
dYH  dxdYH  ~  fadY*  =  R*ldYHTdY}i 


(2.16) 


(2.17) 


This  can  be  reduced  to  an  ordinary  differential  equation  by  making  the  transfor¬ 
mation: 

*  =  FM\[srl  ”  =  W?  (2-18) 

If  we  also  assume  that  viscosity  varies  linearly  with  temperature  we  obtain 

2  Fm  +  FF"  =  0  (2.19) 

where  the  prime  here  denotes  differentiation  with  respect  to  rj.  Boundary  conditions 
are  that  F(0)  =  0,  the  dividing  streamline,  and  the  free  stream  velocities  F'( oo)  =  1 
and  F'(-oo)  =  Ui.  Equation  (2.19)  can  then  be  solved  by  shooting. 

The  shooting  procedure  is  simplified  by  using  the  invariance  of  the  equation  to 
the  transform  G(£  =  art)  =  F(rj)/a ,  reducing  the  shooting  parameters  from  2  to  1. 
The  method  is  as  follows. 

(a)  Guess  F'( 0)  and  F"( 0)  and  shoot  to  +00. 

(b)  Evaluate  the  constant  a  and  transform  the  equation  so  that  G  satisfies  the 
upper  boundary  condition  exactly. 


(c)  Shoot  to  — oo  and  compare  G'{  — oo)  with  U2. 

(d)  Choose  new  F"(0)  and  iterate. 

A  fourth  order  Runge-Kutta  scheme  was  used  for  the  integrations  and  secant 
iteration  was  used  to  vary  F"( 0)  until  G^-oo)  converged  to  U2.  Once  the  mean 
velocity  profile  G'(Yff)  is  known  the  temperature  can  be  found  using  (2.13).  The 
shooting  cannot  proceed  to  infinity,  so  a  cutoff  at  Yfj  =  20  was  used,  after  checking 
that  this  had  negligible  effect  on  the  results.  The  last  step  in  the  procedure  is 
to  reverse  Howarth’s  transform  and  convert  back  to  the  physical  y  coordinate  by 
integrating: 


^=T  V(0)  =  0 


(2.20) 


and  normalizing  the  y  coordinate  by  the  vorticity  thickness,  defined  by 


1-U2 


(2.21) 


Papamoschou  [1986]  noted  that  the  results  presented  by  Gropengiesser  [1970] 
showed  a  larger  than  expected  sensitivity  to  density  ratio.  Both  the  current  calcu¬ 
lations,  and  those  of  Ragab  &  Wu  [1988]  show  that  the  work  of  Gropengiesser  was 
not  in  error.  The  reason  for  the  discrepancy  lies  in  the  normalization  of  the  mean 
profiles.  Gropengiesser  normalized  the  thickness  by  the  momentum  thickness  in  the 
Yj{  domain.  Thus  his  plots  for  different  density  ratios  are  normalized  differently, 
according  to  the  effect  of  density  on  the  transform  Yjj  — ►  y.  To  properly  assess 
the  effect  of  density  ratio  on  amplification  rate  one  has  to  normalize  by  a  consis¬ 
tent  thickness  parameter.  The  vorticity  thickness  was  chosen  since  Monkewitz  and 
Huerre  [1982]  found  that  vorticity  thickness  rather  that  momentum  thickness  gave 
direct  proportionality  between  spatial  amplification  rate  and  mixing  layer  growth 
rate. 

Rather  than  fit  the  profiles  to  a  generalized  hyperbolic  tangent,  as  done  by 
Gropengiesser  [1970],  it  was  decided  to  use  the  computed  profiles  directly.  A  cubic 
spline  was  fitted  through  the  integration  points  and  this  was  used  as  the  base  profile 
for  stability  calculations. 
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2.1.2  Linearized  Disturbance  Equations  and  Shooting  Method 


The  inflectional  instability  of  the  mixing  layer  provides  an  inviscid  instability 
mechanism,  and  the  only  effect  of  viscosity  is  to  damp  the  growing  disturbances 
(Betchov  and  Szewczyk  [1963]).  It  was  therefore  decided  to  solve  the  simpler  inviscid 
stability  problem.  The  starting  equations  are  the  Euler  equations,  obtained  by 
dropping  the  viscous  and  heat  conduction  terms  from  the  Navier-Stokes  equations. 
In  dimensional  form  these  equations  for  continuity,  momentum  and  energy  are  as 
follows  (White  [1974]): 


.«/>*  .a<  _ 

at'  +  “■  ax'  +  f  ax'  0 

(2.22) 

ap- 

P  \dt*  >  dx*)  dx * 

J  * 

(2.23) 

.dh*\  dp *  .dp* 

P  \dt*+  U*  dx* )  dt*  +  “*  dx * 

(2.24) 

Using  the  continuity  equation,  the  perfect  gas  equation  (2.4),  and  the  definitions  of 
h *  and  qf  from  the  previous  subsection,  the  energy  equation  can  be  rewritten  as: 


.(dT*  .dT*\  ..  .ou- 

P  (dt*  +  U'  dx*)~  P  ^  ^dx* 


du* 


(2.25) 


Non-dimensionalization  is  obtained  by  the  same  method  as  in  the  previous  sec¬ 
tion,  equation  (2.5).  The  non-dimensional  equations  for  continuity,  momentum  and 
energy  are: 


dp  dp  du . 

dt  *  di{  ^  dx,- 

(2.26) 

/dut  du,\  1  dp 

\  dt  dij  )  qMj  dx± 

(2.27) 

dT  dT  \  .  .  du{ 

dT  +  u^)  =  -p('y-1fe 

(2.28) 
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and  the  non-dimensional  perfect  gas  equation  is 

p  =  pT  (2.29) 

We  introduce  the  following  decomposition: 

u  =  u  +  u',v  =  v',w  =  w',T  =  T  +  T',p  =  p-\-  p',p  =  1  +  p1  (2.30) 

where  an  overbar  denotes  a  mean  quantity  and  a  prime  a  small  disturbance.  It  is 
assumed  that  the  stability  analysis  may  be  conducted  neglecting  the  slow  streamwise 
variation  of  the  mean  flow,  t.e.  u  =  tZ(y),  T  =  T(y)  etc.  The  mean  velocity,  density 
and  temperature  profiles  need  to  be  specified.  The  mean  pressure  is  constant  and 
non-dimensionally  unity. 

Under  the  parallel-mean-flow  assumption,  the  linearized  disturbance  equations 
have  coefficients  which  are  independent  of  x ,  z  and  t.  Hence  the  solutions  are  expo¬ 
nentials  in  these  independent  variables  and  disturbances  have  the  form  of  travelling 
waves 

(u',  v\  w',  T\  p\  p')  =  (u,  t>,  u>,  f,  p,  p)e'{<*x+0*-ut)  (2.3i) 

where  u,  v,  w,  T,  p,  p  are  complex  eigenfunctions  depending  only  on  the  y  coordinate. 
In  equation  (2.31)  u  is  the  frequency  and  a  and  (3  are  wavenumbers  in  the  streamwise 
(x)  and  spanwise  (z)  directions  respectively.  The  angle  6  of  a  disturbance  is  given 

by 

tan  0  —  (3/ar  (2.32) 

where  ar  is  the  real  part  of  a.  The  wavenumber  /?  has  to  be  real,  since  we  require 
disturbances  not  to  amplify  for  z  — »  ±oo.  The  form  of  a  and  u;  depends  on  the 
particular  problem.  For  the  temporal  stability  problem  disturbances  grow  in  time 
and  not  in  space,  so  a  has  to  be  real  and  w  complex.  In  the  spatial  problem  u 
is  real  and  a  complex  so  that  disturbances  grow  in  space,  but  not  in  time.  The 
amplification  rate  of  a  disturbance  is  given  by  u/t-  in  the  temporal  case  and  —a,  in 
the  spatial  case. 

The  pressure  perturbation  can  be  easily  found  from  the  density  and  temperature 
using  the  linearized  form  of  the  perfect  gas  law  (2.29): 

p'  =  pT'  +  p'T  (2.33) 
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After  substituting  equations  (2.30)  and  (2.31)  into  the  governing  equations  (2.26)- 
(2.29)  and  linearizing,  we  obtain  the  following  linearized  equations  for  continuity, 
three  components  of  momentum,  energy  and  the  perfect  gas  law: 


pt(au  —  w)  +  vDp  +  p[»(au  +  fiw)  +  Dv]  =  0 


(2.34) 


p[t(au  -  w)u  +  vDu]  =  ~ 


(2.35) 


pi(au  —  u)v  = 


-Dp 
1 M* 


(2.36) 


pt(au  —  a;)tD  = 


-ifip 


(2.37) 


p[t(au  —  w)T  +  vDT]  =  —(7  —  l)[*(au  +  pw)  +  Dv] 


(2.38) 


p  =  pT  +  pT 


(2.39) 


where  D  represents  the  operator  d/dy . 

These  equations  can  be  reduced  to  a  set  of  two  equations  as  follows.  First  we 
multiply  equation  (2.34)  by  T  and  then  add  it  to  (2.38),  using  the  fact  that  differ¬ 
entiating  (2.29)  with  respect  to  y  gives  TDp  +  pDT  =  0.  This  gives 


*(au  —  u)p  =  —7(1 (au  +  0w)  +  Dv] 

Now  u  and  u>  can  be  eliminated  using  equations  (2.35)  and  (2.37),  giving 


(2-40) 


(au  —  u)Dv  —  avDu  = 


•A  2 

ipaLg 


(2.41a) 


where 


a2  +  /?2  A/?(au  — w)i 


Equations  (2.36)  and  (2.41)  form  a  set  of  two  equations  for  t;  and  p.  These  can 
be  further  reduced  by  following  the  method  of  Gropengiesser  [1970],  who  defined  a 
new  variable  x: 

tap 

The  equation  to  be  integrated  for  x  then  becomes: 


dx  _  a2(u-^/a)  _  x(xg  +  77u) 
dy  T  (u  -  w/a) 


(2.43) 


The  boundary  conditions  are  obtained  by  considering  the  asymptotic  behavior  of 
disturbances  in  the  freestreams.  In  the  freestream  Du  is  zero  and  equations  (2.36) 
and  (2.41)  can  be  written  as 

D2p  =  pa2gp  (2.44a) 

D2  v  =  pa2gv  (2.446) 

The  general  solutions  to  these  equations,  vanishing  for  y  — *  ±oo  are 

p  =  aie*™  (2.45a) 

v  =  a2eTqy  (2.456) 

where  a*  and  a2  are  constants  and  q  is  ~pa2g.  Note  that  p  and  v  decay  with  the 
same  complex  exponential  in  the  freestream.  Therefore,  the  new  variable  x,  formed 
as  a  ratio  of  these  solutions,  must  go  to  a  constant  as  y  — >  ±oo.  Setting  dx/dy  =  0 
and  Du  =  0  in  equation  (2.43)  leaves 


X(y  =  ±oo)  =  ->  — t=M  2.46 

VffT 

The  numerical  solution  procedure  is  iterative.  First,  a  guess  is  made  of  the 
eigenvalue.  For  spatial  analysis  u>  is  specified  and  a  guessed,  whereas  for  temporal 
analysis  a  is  specified  and  u>  guessed.  Knowing  the  eigenvalue,  we  can  evaluate  x  in 
the  freestreams  using  equation  (2.46).  Then  we  integrate  equation  (2.43)  from  each 
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of  the  freestreams  to  the  centerline  of  the  mixing  layer  at  y  =  0.  At  the  centerline  the 
value  of  x  computed  by  integrating  from  the  upper  freestream,  x+(0),  is  compared 
with  the  value  computed  by  integrating  from  the  lower  freestream,  x_(0)-  Then  a 
new  eigenvalue  is  chosen,  and  iterated  upon  until  the  eigenvalue  has  converged  to  a 
specified  tolerance.  The  shooting  method  was  implemented  using  subroutines  from 
Press  ct  al.  [1986].  The  integrations  were  carried  out  using  a  variable  step  fifth 
order  Runge-Kutta  scheme,  and  the  iteration  was  achieved  by  a  Newton-Raphson 
method.  An  error  control  of  10~6  was  used  for  the  integrations,  and  the  iteration 
continued  until  the  eigenvalue  converged  to  10~7.  Single  precision  arithmetic  was 
sufficient  for  most  computations.  However,  double  precision  was  required  for  some 
of  the  weakly  amplified  supersonic  modes  of  instability  (section  2.2.4). 

Generally  the  integrations  were  performed  with  y  as  the  independent  variable, 
starting  the  integrations  10  vorticity  thicknesses  away  from  the  centerline.  Instead 
of  solving  the  problem  on  the  domain  [— oo,  oo]  it  is  possible  (Gropengiesser  [1970]) 
to  reduce  the  domain  of  integration  to  [C/2, 1]  by  transforming  the  independent 
variable  in  equation  (2.43)  from  y  to  u.  This  method  was  also  implemented  and 
was  found  to  be  slightly  quicker  than  the  method  integrating  in  y,  and  gave  exactly 
the  same  results.  However,  there  were  problems  with  this  method  for  wake  flows  and 
for  mixing  layer  profiles  where  the  base  flow  was  not  an  analytic  function.  Results 
from  compressible  wake  calculations  using  the  above  methods  can  be  found  in  Chen 
et  al.  [1989].  The  procedure  for  the  wake  calculations  was  similar.  However,  for  the 
symmetric  wake  mode  v  at  the  centerline  is  zero,  so  x  defined  by  equation  (2.42) 
goes  to  infinity.  This  was  remedied  by  working  with  a  new  variable  1/x- 

Once  the  eigenvalues  have  been  found  the  eigenfunctions  are  calculated  by  in¬ 
tegrating  equations  (2.36)  and  (2.41)  out  from  the  centerline  into  the  freestream. 
Initial  conditions  are  calculated  from  the  solution  for  x(0),  choosing  the  phase  of 
the  eigenfunctions  so  that  v(0)  =  14-  *0.  From  (2.42)  this  means  that  p(0)  is  given 
by: 

pr(0)  =  iMlxM  (2.47a) 

P,(0)  =  0)  (2.476) 

After  computation  the  eigenfunctions  are  renormalized  without  phase  change  so 
that  the  magnitude  of  u  is  1.  In  all  cases  this  renders  the  other  components  of  the 
eigenfunction  less  than  1. 
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2.1.3  Validation 


The  first  method  of  validation  was  to  compare  results  obtained  with  the  shooting 
method  described  in  the  previous  section  with  results  from  a  code  for  direct  solution 
of  the  linearized  equations,  developed  in  collaboration  with  J.  H.  Chen  and  described 
in  Appendix  A.  This  is  an  ideal  situation  for  checking  numerical  methods.  When  two 
completely  different  methods  give  the  same  results  one  is  confident  that  the  solution 
is  correct.  The  superior  performance  of  the  shooting  method  for  weakly  amplified 
disturbances  led  to  its  being  used  for  all  the  stability  calculations  presented  here. 

In  addition,  results  were  checked  against  published  data  including  the  graphs  of 
Gropengiesser  [1970],  and  at  low  Mach  number  against  results  for  incompressible 
flow  by  Michalke  [1965a, b],  Monkewitz  and  Huerre  [1982]  and  Lowery  and  Reynolds 
[1986].  Comparisons  are  presented  below  in  tables  2.1  and  2.2. 

Table  2.1  Comparison  of  temporal  results  at  M\  =  0.01  with  Michalke  [1965] 


[ 

I 

^  Table  2.2  Comparison  of  spatial  results  at  =  0.02 

■  with  Lowery  and  Reynolds  [1986] 


I 
l 

|  2.2  Results 

L  Results  are  presented  for  a  variety  of  mixing  layers.  In  non-dimensional  terms 

I  (equation  (2.5))  these  flows  are  characterized  by  U2,  the  velocity  ratio  and  P2,  the 

density  ratio.  Since  pressure  is  assumed  uniform  we  can  also  write  P2  =  I/T2. 


ur 

a  Lowery 

a  current 

■a | 

0.88869,-0.12850* 

0.88891,-0.12850* 

0.43110,-0.09913* 

0.43110,-0.09913* 

Hi 

0.20908,  -0.05860* 

0.20908,  -0.05860* 

ar 

Michalke 

u>{  current 

0.2 

0.06975 

0.06974 

0.4 

0.09410 

0.09409 

0.6 

0.08650 

0.08649 

0.8 

0.05388 

0.05386 

2.2.1  Low  Mach  Number  Results 

It  was  shown  by  Monkewitz  and  Huerre  [1982]  that,  when  spatial  theory  and  the 
Blasius  mixing  layer  profile  are  used,  the  maximum  amplification  rate  |o^|max  is 
proportional  to  A  =  ( U ^  —  £/£ )/ (^i  +  ^so>  the  best  fit  through  the  experi¬ 
mental  data  for  the  incompressible  mixing  layer  with  uniform  density,  compiled  by 
Brown  and  Roshko  [1974],  is  a  straight  line  dS/dx  proportional  to  A.  Thus,  for  the 
incompressible  mixing  layer  with  constant  density  it  appears  that  |c^|max  from  the 
linear  theory  is  proportional  to  dSjdx  from  experiments  (Morkovin  [1988]). 

In  this  section  we  test  the  postulate  that  the  relation  |o^[max  ~  dS  jdx  applies  to 
all  mixing  layers,  with  any  velocity  ratio,  density  ratio  and  Mach  number,  with  the 
following  provisos: 

(a)  use  spatial  stability  theory  since  this  is  a  convectively  unstable  flow  (Huerre 
and  Monkewitz  [1985]) 

(b)  use  a  solution  of  the  laminar  boundary-layer  equations  as  the  base  flow 

(c)  normalize  the  profiles  by  the  vorticity  thickness. 

The  effect  of  the  mean  velocity  profile  can  be  important.  A  comparison  between 
hyperbolic  tangent  velocity  profiles  and  the  boundary-layer  solution  is  shown  for 
two  different  velocity  ratios  in  figure  2.1.  In  each  case  the  densities  of  the  free- 
streams  are  equal  (/>2  =  1)-  The  profiles  are  all  normalized  with  their  vorticity 
thickness  and  shifted  so  that  u  =  (l  +  C^2)/2  at  y  =  0.  The  most  amplified  spatial 
instability  wave  was  computed  from  these  profiles.  The  percentage  difference  in 
amplification  rate  of  the  case  with  a  hyperbolic  tangent  profile  relative  to  the  case 
with  a  boundary-layer  profile  is  shown  on  table  2.3. 

Table  2.3  Percentage  difference  in  lo^lmax  between  tanh  and 
laminar  velocity  profiles  at  Mi  =  0.1  and  P2  =  1 
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The  largest  effect  was  at  U2  =  0  where  the  wave  with  the  highest  amplification 
rate  was  16%  more  rapidly  amplified  on  the  hyperbolic  tangent  mean  flow  than  on 
the  boundary- layer  mean  flow.  At  higher  values  of  U2  the  effect  was  much  smaller. 

To  investigate  the  effect  of  profile  shape  some  calculations  were  made  with  a  pro¬ 
file  constructed  from  a  hyperbolic  tangent  on  the  low  speed  side  and  the  boundary- 
layer  solution  on  the  high  speed  side.  Such  profiles  are  smooth  to  the  eye  but 
derivatives  are  discontinuous  at  the  centerline.  No  numerical  problems  were  en- 
countrered  but  results  have  to  be  treated  with  caution.  It  was  found  that  a  ‘fuller’ 
velocity  profile  on  the  high  speed  side  was  stabilizing,  and  from  the  reverse  cal¬ 
culation  that  a  longer  tail  on  the  low  speed  side  was  also  stabilizing.  The  latter 
effect  was  larger  by  a  factor  of  about  3.  This  kind  of  argument  may  explain  why 
tripping  the  high  speed  boundary-layer,  making  it  turbulent  and  hence  having  a 
much  ‘fuller’  mean  velocity  profile  reduces  the  growth  rate  of  the  developed  mixing 
layer  by  approximately  30%  (Browand  and  Latigo  [1979],  Mungal  et  al.  [1985]). 

The  differences  in  the  velocity  profile  are  generally  small  when  the  two  free- 
streams  have  equal  densities,  but  can  become  very  large  when  there  is  a  large 
density  ratio.  The  effect  of  density  ratio  on  the  mean  velocity  and  density  profiles 
at  U2  =  0.5  and  p2  =  1/7, 1,7  is  shown  in  figure  2.2.  The  profiles  for  P2  =  1/7 
and  P2  =  7  are  very  different  from  a  hyperbolic  tangent.  Clearly  in  this  situation  a 
hyperbolic  tangent  would  be  a  poor  choice  for  the  mean  velocity  profile.  The  most 
unstable  case  corresponds  to  P2  =  7,  which  is  the  least  ‘full’  profile  on  the  high 
speed  side  and  the  shortest  tail  on  the  low  speed  side,  in  agreement  with  the  above 
arguments. 

The  effect  of  density  and  velocity  ratios  on  the  amplification  rate  of  the  most 
amplified  spatial  instability  wave  (|a,  |max)  is  shown  in  figure  2.3a  for  the  boundary- 
layer  mean  velocity  profile  and  on  figure  2.36  for  the  hyperbolic  tangent  mean 
velocity  profile  The  growth  rate  is  plotted  against  A  for  three  different  density 
ratios,  as  done  by  Bogdanoff  [1984]  and  Dimotakis  [1986].  Experimental  points  from 
Brown  and  Roshko  [1974]  are  also  plotted,  both  the  original  vorticity  thicknesses  as 
well  as  density  thicknesses.  The  density  thicknesses  were  evaluated  by  Bogdanoff 
[1984]  by  joining  the  20%  and  80%  points  on  the  density  profile  and  measuring  the 
distance  between  the  points  where  this  line  intercepts  the  free-stream  density.  The 
agreement  between  the  linear  amplification  rate  and  the  experimental  growth  rate 
is  remarkable,  especially  when  mixing  layer  thicknesses  based  on  the  experimental 
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mean  density  profile  are  used.  For  these  and  other  experiments  the  ratio  61  /  |a,| 
where  6'  =  d6/dx,  is  shown  on  table  2.4  below. 


Table  2.4  Relationship  between  6'  from  experiments  and  |aj 

from  linear  theory 


Experiment 

u2 

P2 

Brown  k  Roshko  1974 

Brown  k  Roshko  1974 

Brown  k  Roshko  1974 
Fiedler  [1974] 

Bogdanoff  [1984] 

Dimotakis  k  Brown  [1976] 

i/Vi 

l/v/7 

1/7 

0.0 

0.0 

0.19 

7 

1/7 

7 

1.09 

0.2 

1.0 

From  these  results,  we  conclude  the  following  relations  between  linear  amplifica¬ 
tion  and  experimental  growth  rate  (±20%): 


=  O^KIr 


(2.48a) 


=  0.4510,1, 


(2.486) 


Comparison  between  figures  2.3a  and  2.36  shows  that  the  linear  results  based 
on  a  hyperbolic  tangent  mean  velocity  profile  do  not  show  the  correct  trends.  At 
A  =  0.3,  for  example,  the  effect  of  increasing  density  ratio  with  the  hyperbolic 
tangent  mean  velocity  profile  is  that  first  the  amplification  rate  rises  and  then  it 
falls.  This  can  be  compared  with  the  results  for  the  boundary-layer  mean  flow, 
which  show  a  continuous  rise  in  amplification  rate  with  increasing  density  ratio. 

The  variation  in  |a,|max  over  two  orders  of  magnitude  change  in  temperature 
ratio  72  (recall  T2  =  1  / P2)  is  shown  in  figure  2.4  for  velocity  ratios  f/2  =  0  and 
0.5.  The  case  U2  —  0  was  computed  by  Maslowe  and  Kelly  [1971]  for  a  hyperbolic 
tangent  mean  velocity  and  a  specified,  non  Crocco-Busemann,  mean  temperature 
profile.  They  found  a  peak  amplification  rate,  at  a  density  ratio  of  p2  =  33,  that 
was  64%  higher  than  the  amplification  rate  for  equal  densities.  The  current  results 
for  that  velocity  ratio  show  a  much  stronger  effect  of  density  ratio  on  amplification 
rate,  in  better  agreement  with  the  experiments  of  Davey  and  Roshko  [1972].  At  a 


density  ratio  p2  =  5,  the  amplification  rate  is  three  times  that  of  the  equal-density 
case.  When  U2  =  0  the  amplification  rate  appears  to  become  very  large  when  the 
low  speed  stream  has  a  large  density  compared  to  the  high  speed  stream.  However, 
it  should  be  noted  that  for  the  case  U2  —  0  the  v  component  of  velocity  on  the  low 
speed  side  of  the  layer  will  be  large  relative  to  the  u  component  and  the  boundary- 
layer  assumption  is  no  longer  strictly  valid. 

We  have  seen  that  there  is  a  good  agreement  between  the  linear  theory  and  the 
existing  experimental  data  for  the  effect  of  density  and  velocity  ratio  on  the  mixing 
layer  growth  rate  at  low  Mach  numbers.  It  therefore  appears  reasonable  to  use 
the  linear  theory  over  the  full  range  of  possible  conditions  (not  just  at  experimental 
points)  to  compare  with  some  of  the  models  that  have  been  proposed  for  the  growth 
rate  of  the  mixing  layer  a s  a  function  of  density  and  velocity  ratios.  In  particular, 
two  models  are  compared.  The  first  is  a  form  originally  proposed  by  Brown  [1974], 
and  used  by  Papamoschou  [1986]  to  compute  the  growth  rates  of  incompressible 
mixing  layers  with  the  same  velocity  and  density  ratio  as  his  compressible  mixing 
layer  experiments: 


d6_  ~  (1  -  +  y/Pi) 

dx  (1  +  U2y/P2) 


(2.49) 


The  second  model  is  a  modified  form  of  the  above,  and  was  proposed  from  geomet¬ 
rical  arguments  by  Dimotakis  [1986] 


dx  {l  +  U2Jn)\  1  + 2.9/A  / 


(2.50) 


These  models  will  be  referred  to  as  model  1  and  model  2.  For  comparison  purposes 
the  constant  of  proportionality  is  set  so  that  both  models  coincide  with  the  linear 
growth  rate  prediction  (equation  (2.48))  for  the  equal-density  case.  The  comparison 
between  the  models  and  the  linear  theory  is  shown  in  figure  2.4  for  velocity  ratios 
0  and  0.5.  Agreement  is  generally  not  good.  We  can  also  check  the  models  by 
comparing  the  growth  rate  plotted  against  the  right  hand  sides  of  equations  (2.49) 
and  (2.50)  above.  A  straight  line  would  indicate  agreement  between  the  model  for 
growth  rate  and  the  linear  theory  prediction.  Figures  2.5a  and  2.56  show  this  plot 
for  each  model.  Model  2  is  closer  to  a  straight  line. 
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Both  model  1  and  model  2,  and  the  linear  theory  prediction  of  equation  (2.48) 
agree  well  with  the  available,  but  limited,  experimental  data.  However  the  different 
predictions  of  these  three  models  for  other  velocity  and  density  ratios  suggests  that 
this  is  not  a  cut-and-dried  issue.  In  particular,  the  method  of  normalization  used 
by  Papamoschou  and  Roshko  [1988]  (dividing  the  compressible  mixing  layer  growth 
by  the  (model  1)  prediction  of  growth  rate  of  an  incompressible  mixing  layer  with 
the  same  velocity  and  density  ratio)  has  the  built-in  assumption  that  model  1  is 
correct. 

2.2.2  Oblique  Waves  at  High  Mach  Number 

The  basic  effect  of  compressibility  is  first  considered  for  the  temporal  stability  of 
the  time-developing  mixing  layer  with  equal  free-stream  densities  and  temperatures, 
and  with  a  simple  velocity  profile: 


u  =  tanh(2y) 


(2.51) 


Equation  (1.3)  for  the  convective  Mach  number  can  be  non-dimensionalized  as  fol¬ 
lows: 

1-^2  _  Miy/P2(l  -  U2) 


Mc  = 


Cl  +  C2 


1  +  y/P2 


(2.52) 


In  the  time-developing  reference  frame  U2  =  —  1  and  for  P2  =  1  we  have  M\  =  Mc, 
so  the  convective  Mach  number  is  the  same  as  the  Mach  number  of  each  of  the 
free-streams. 

Figure  2.6a  shows  the  effect  of  increasing  Mach  number  on  the  amplification  rate 
of  two-dimensional  waves.  The  observed  trend  is  the  same  as  found  by  Gropengiesser 
[1970]  for  the  spatial  theory,  including  the  appearance  of  a  second  (supersonic)  mode 
of  instability.  This  mode  is  characterized  by  a  different  phase  speed  cr  =  wr/ar, 
as  evident  from  the  plot  of  wr  shown  in  figure  2.66.  The  effect  of  increasing  Mach 
number  is  to  reduce  the  amplification  rate  and  ultimately  stabilize  the  subsonic  two- 
dimensional  mode.  Only  the  emergence  of  the  second  mode  keeps  the  mixing  layer 
unstable  to  two-dimensional  disturbances  at  high  Mach  numbers.  The  subsonic 
(first)  mode  is  stationary  in  the  time-developing  mixing  layer  reference  frame,  while 
the  second  mode  travels  to  the  left  or  right  (see  section  2.2.4). 
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The  amplification  rates  of  oblique  waves  at  angles  of  0°,  30°,  and  60°  are  shown 
in  figures  2.7  through  2.9  for  Mach  numbers  0.01,  0.8,  and  1.2  respectively.  At 
Mi  =  0.01  the  most  amplified  wave  is  the  two-dimensional  wave  ( 6  =  0°),  as 
expected  from  Squires  theorem,  which  applies  only  to  incompressible  flow  (Squire 
[1933]).  At  Mi  =  0.8  it  is  found  that  waves  of  0°,  30°  and  60°  are  all  about  equally 
amplified,  and  by  Mi  =  1.2  the  wave  at  60°  is  substantially  more  amplified  than 
the  less  oblique  waves. 

The  increasing  obliquity  of  the  most  amplified  waves  is  better  illustrated  in  fig¬ 
ures  2.10  and  2.11.  In  these  figures  the  amplification  rate  is  plotted  against  9  for 
various  Mach  numbers,  where  for  each  Mach  number  the  wavelength  is  fixed  at 
the  most  amplified  wavelength  (including  oblique  waves).  Figure  2.10  is  the  plot 
for  the  time-developing  mixing  layer  from  above,  whereas  figure  2.11  is  for  a  spa¬ 
tially  developing  mixing  layer  with  —  1.0  and  =  0.5,  and  with  a  compressible 
laminar  boundary-layer  solution  as  the  base  flow.  The  plots  are  very  similar  for 
the  two  cases,  indicating  that  for  mixing  layers  with  equal  densities  the  hyperbolic 
tangent  velocity  profile  can  be  reliably  used  to  compute  the  fundamental  effects 
of  compressibility.  In  each  case  the  curves  split  into  two  regimes.  For  Mc  <  0.6 
the  two-dimensional  wave  is  always  the  most  amplified  wave.  Above  Mc  =  0.6  a 
three-dimensional  wave  of  increasing  obliquity  is  most  amplified.  The  second  mode 
is  amplified  at  high  Mach  number,  but  never  more  than  the  oblique  first  mode.  The 
cusp  in  the  plots  corresponds  to  the  transition  as  one  mode  becomes  more  amplified 
than  the  other. 

The  angle  of  the  most  amplified  distu.  bance  was  determined  empirically  to  satisfy 

Mc  cos  0  as  0.6  (2.53) 

This  means  that  the  Mach  number  perpendicular  to  the  wave  crest  is  approximately 
0.6.  This  might  be  considered  similar  to  the  case  of  swept-back  airfoils,  where  the 
key  Mach  number  is  Moo  cos  9.  In  the  mixing  layer,  waves  at  0  =  90°  are  not 
amplified  and  two-dimensional  waves  ( 0  =  0°)  have  growth  rates  that  are  strongly 
reduced  by  compressibility  effects  (see  section  2.2.4  and  Chapter  4).  The  most 
amplified  wave  has  to  be  somewhere  in  between. 

Table  2.5  shows  values  of  Mc  cos0  for  spatially-developing  mixing  layers  con¬ 
structed  in  three  different  ways: 


29 


(a)  Equal  temperatures  and  fixed  velocity  ratio:  T2  =  1  and  U2  =  0.5 

(b)  Equal  stagnation  enthalpies  (dimensionally  we  have  H*  =  c*T*  +  u*2/2,  non- 
dimensionally  H  —  T / ((7  —  1  )M2 )  -f  u2/2  )  and  zero  velocity  on  the  low-speed 
side:  H2  =  1  and  U2  =  0.  This  case  corresponds  to  the  earliest  experiments 
on  the  compressible  mixing  layer. 

(c)  Fixed  Mach  numbers  Mi  =  2.0  and  M2  =  1.0.  In  this  case  the  convective 
Mach  number  is  varied  by  changing  the  ratio  of  stagnation  enthalpies  and  the 
velocity  ratio. 

In  each  case  a  boundary-layer  mean  flow  was  used,  for  comparison  with  experi¬ 
ments.  The  value  of  MccosO  is  consistently  between  0.58  and  0.59  for  case  (a)  but 
shows  some  deviation  at  high  Mach  number  for  the  cases  (b)  and  (c). 


Table  2.5  Variation  of  Mc  cos  0  for  oblique  waves 


T2  = 

1,  u2  =  0.5 

H2  =  1,U2  =  0 

M\  =  2,  M2  =  1 

Mc 

Mc  cos  0 

Me 

Me  cos  0 

Mc 

Mc  cos  0 

0.6 

0.587 

0.571 

0.599 

0.8 

0.585 

1.122 

0.623 

Era 

0.601 

1.0 

0.588 

1.311 

0.734 

K  jm 

0.628 

0.582 

BBSs 

0.651 

Mm 

0.581 

1.045 

0.699 

wBm 

0.586 

1.107 

0.769 

2.2.3  Convective  Mach  Number 

The  variation  of  the  amplification  rate  of  the  most  unstable  mode  with  Mach 
number  is  shown  in  figure  2.12a  for  the  temporal  stability  of  the  time-developing 
mixing  layer.  The  curve  of  the  most  amplified  two-dimensional  wave  is  found  by 
varying  the  wavenumber  ar  until  a  maximum  is  found,  keeping  0  fixed  at  0°.  The 
curve  for  oblique  waves  is  found  by  varying  both  wavenumber  and  angle  until  a 
maximum  is  obtained.  Above  M\  =  0.6  the  three-dimensional  waves  are  the  most 
rapidly  amplified  and  the  curve  of  the  most  amplified  oblique  disturbance  is  a  much 
better  fit  to  the  existing  experimental  data  (Papamoschou  and  Roshko  [1988],  see 
also  figure  2.14)  than  the  two-dimensional  curve.  Figure  2.126  shows  the  extension 
of  the  curve  for  oblique  waves  to  very  high  Mach  numbers.  It  can  be  observed  that, 
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at  least  for  the  case  of  equal  densities,  the  growth  rate  continues  to  decrease  as 
Mach  number  is  increased. 

To  check  the  convective  Mach  number  concept,  and  to  provide  data  for  compar¬ 
ison  with  experiment,  curves  of  the  maximum  spatial  amplification  rate,  |o,  jmax, 
against  convective  Mach  number  were  compiled.  In  each  case  the  compressible 
boundary-layer  equations  were  solved  for  the  mean  flow,  and  the  most  amplified 
disturbances  (as  a  function  of  both  frequency  and  angle)  were  determined.  These 
were  normalized  with  the  amplification  rate  of  the  most  amplified  disturbance  at 
Mi  =  0.1,  with  the  same  velocity  and  density  ratio.  Figure  2.13  shows  the  graphs 
of  peak  amplification  rate  against  Mach  number,  using  the  three  different  methods 
of  varying  the  convective  Mach  number  described  in  the  previous  section.  A  good 
collapse  of  the  data  with  Mc  is  obtained  for  Mc  <  0.8,  but  there  is  some  divergence 
at  high  convective  Mach  numbers. 

It  appears  from  recent  work  by  Papamoschou  [1989]  that  the  fundamental  idea 
behind  the  convective  Mach  number  concept,  of  a  large-scale  structure  convect- 
ing  in  the  flow  and  ‘seeing’  the  relative  Mach  number  of  the  free-streams,  may  be 
flawed,  since  actual  measured  convective  velocities  do  not  agree  with  the  theory. 
An  alternative  viewpoint  is  that  Mc,  as  defined  by  equation  (1.3),  (i.e.  the  velocity 
difference  divided  by  twice  the  average  sound  speed)  is  just  a  dimensionless  param¬ 
eter  that  can  include  only  the  first-order  effects  of  compressibility.  With  this  latter 
viewpoint  we  do  not  expect  perfect  collapse  of  growth  rates  at  high  Mach  numbers. 

It  was  mentioned  in  section  2.2.1  that  a  source  of  error  in  the  normalization  used 
by  Papamoschou  and  Roshko  [1988]  is  the  model  for  incompressible  growth  rate, 
equation  (2.49).  An  alternative  method  it  to  normalize  experimental  data  by  the 
growth  rate  of  the  most  amplified  spatial  instability  wave,  |at|max,  at  low  Mach 
number.  This  is  shown  on  figure  2.14  for  the  available  pitot  thickness  data  from  Pa¬ 
pamoschou  [1986],  together  with  a  typical  curve  for  the  linear  theory  with  T 2  =  1.0. 
The  linear  theory  curve  has  been  anchored  to  be  0.6  at  M\  =  0,  using  equation 
(2.48a)  and  assuming  that  density  thickness  and  pitot  thickness  are  equal.  Unfor¬ 
tunately  the  large  spread  in  the  experimental  data  prevents  a  definitive  conclusion 
on  the  performance  of  this  method  of  normalization. 
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2.2.4  Supersonic  Instability  Modes 

We  have  seen  that  the  two-dimensional  instability  mode  that  is  most  amplified 
at  low  Mach  number  becomes  stable  at  high  Mach  number.  Previous  researchers 
(Lessen  tt  al.  [1965,1966],  Gropengiesser  [1970])  have  found  other  modes  of  instabil¬ 
ity  at  high  Mach  number  which  keep  the  mixing  layer  unstable  in  two-dimensions  at 
any  Mach  number.  These  modes  have  been  described  as  ‘radiating  vorticity  modes’ 
by  Mack  [1989],  and  are  distinct  from  the  ‘Mack  modes’  found  in  supersonic  wall 
boundary  layers,  confined  shear  layers,  and  compressible  wakes,  which  require  the 
presence  of  a  region  of  trapped  subsonic  flow  relative  to  the  free-stream.  The  ra¬ 
diating  vorticity  modes  Eire  supersonic  with  respect  to  one  of  the  free-streams  and 
radiate  energy  into  that  stream.  Two  such  modes  exist  in  the  mixing  layer,  one 
mode  supersonic  relative  to  the  low-speed  stream,  and  the  other  mode  supersonic 
relative  to  the  high-speed  stream. 

Jackson  and  Grosch  [1988]  investigated  the  supersonic  modes  for  the  mixing 
layer  described  by  a  hyperbolic  tangent  profile  in  Howarth- transformed  space  (see 
equation  (2.16)).  In  the  present  work  a  similar  case  to  Jackson  and  Grosch  is 
studied,  but  using  a  solution  to  the  boundary-layer  equations  as  the  base  flow.  The 
flow  chosen  has  U2  —  0  and  T2  =  1  and  spatial  stability  theory  is  used.  Figure  2.15 
shows  the  variation  in  growth  rate  of  the  most  amplified  two-dimensional  modes 
with  Mach  number,  in  the  interesting  region  around  a  free-stream  Mach  number 
Mi  =  2  (i.e.  around  Mc  =  1  since  for  T2  =  1  and  U2  =  0  the  convective  Mach 
number  is  always  one  half  the  Mach  number  of  the  high  speed  stream).  Figure  2.16 
shows  the  phase  speeds  of  these  modes,  which  led  Jackson  and  Grosch  to  classify 
them  as  a  fast  supersonic  mode,  supersonic  relative  to  the  low  speed  stream,  and 
a  slow  supersonic  mode,  supersonic  relative  to  the  high  speed  stream.  The  various 
modes  are  plotted  as  a  function  of  frequency  in  figures  2.17  and  2.18  at  Mi  =  2.2. 
Figure  2.17  shows  the  growth  rate  and  figure  2.18  the  phase  speed.  From  these 
plots  we  see  that  the  most  amplified  mode  is  the  fast  mode.  The  second  peak  at 
the  left  of  the  curve  for  the  fast  mode  is  the  remnant  of  the  subsonic  mode,  which 
was  the  dominant  mode  at  low  Mach  number.  As  Mach  number  is  increased  this 
mode  becomes  less  and  less  amplified,  and  the  peak  moves  to  longer  and  longer 
wavelengths.  From  vortex  sheet  instability  theory  (Blumen  et  al.  [1975])  we  expect 
this  mode  to  finally  become  stable  at  Mi  =  2\/2. 
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The  radiative  nature  of  these  instability  waves  becomes  apparent  when  the  eigen¬ 
functions  for  pressure  are  plotted.  The  same  case  as  above  is  chosen,  with  M\  —  2.2. 
Figure  2.19a  shows  the  pressure  eigenfunction  for  the  fast  supersonic  mode,  radiat¬ 
ing  into  the  low  speed  side,  and  figure  2.196  shows  the  slow  supersonic  mode,  which 
radiates  into  the  high-speed  side.  Figure  2.19c  shows  the  subsonic  non-radiating 
mode.  The  eigenfunctions  shown  in  figure  2.19  have  the  form  of  damped  waves, 
since  we  are  considering  amplified  waves,  and  the  disturbances  away  from  the  cen¬ 
terline  were  created  at  an  earlier  time  and  then  propagated  into  the  free-stream. 
Jackson  and  Grosch  [1988]  only  showed  eigenfunctions  for  the  neutral  instability 
modes,  which  do  not  decay  in  the  free-stream. 

Which  mode  is  most  amplified  depends  upon  the  temperature  ratio.  At  a  tem¬ 
perature  ratio  of  2  the  fast  mode  is  the  dominant  mode,  while  at  a  temperature 
ratio  of  0.5  the  slow  mode  is  dominant.  Unlike  the  subsonic  mode,  increasing  the 
angle  of  the  supersonic  mode  disturbance  does  not  increase  its  amplification  rate, 
and  for  highly  oblique  disturbances  these  modes  are  stable. 

The  supersonic  modes  are  very  interesting  from  a  physical  perspective.  However, 
it  should  be  remembered  that  the  most  amplified  waves  in  the  flow  are  the  oblique 
modes  of  the  subsonic  instability,  and  since  these  are  amplified  more  rapidly  by 
a  factor  of  3  or  4  we  expect  the  resulting  flow  to  be  dominated  by  the  oblique 
waves.  The  supersonic  modes  will  come  into  play  only  if  there  is  very  strong  two- 
dimensional  forcing  of  the  mixing  layer. 

2.2.5  Eigensolution  Structure 

Eigenfunctions  from  the  linear  theory  are  used  in  following  chapters  as  inputs 
to  direct  numerical  simulations.  However  they  are  interesting  in  their  own  right, 
providing  important  clues  to  the  structure  and  physics  of  the  large  scale  motions 
which  develop  from  the  linear  instability.  In  this  section  eigenfunctions  from  two- 
dimensional  stability  calculations  for  the  time-developing  mixing  layer  are  presented 
to  illustrate  the  effects  on  compressibility  on  the  linear  eigenfunctions. 

Figures  2.20  and  2.21  show  the  eigenfunctions  u,t >,p,T  as  functions  of  the  y 
coordinate  for  Mi  =  0.01  and  M\  =  0.6  respectively.  The  main  effect  is  the  rise 
in  importance  of  the  density  and  temperature  eigenfunctions.  These  keep  the  same 
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basic  shape  at  Mi  =  0.6  but  are  increased  by  four  orders  of  magnitude  relative  to 
the  case  M\  =  0.01. 

Additional  information  can  be  found  by  generating  contour  plots  of  the  eigen¬ 
functions.  For  example  the  real  part  of  the  u  velocity 

u  =  u  -f  A  Real[uetaz)  (2.54) 

can  be  found  over  one  wavelength  in  the  x  direction,  and  similarly  derivatives  of 
the  flow  variable  can  be  easily  found.  Results  are  shown  on  figures  2.22a-g  for 
various  flow  quantities  at  M\  =  0.6.  Note  that  the  y  axis  has  been  stretched  to 
better  illustrate  the  structure.  The  amplitude  of  the  disturbance  A  is  chosen  to  be 
0.5  for  the  spanwise  vorticity  u>2,  as  well  as  for  u >2/p  and  p,  to  better  illustrate  the 
structure,  which  would  otherwise  be  dominated  by  the  mean  flow.  The  remaining 
plots  use  a  disturbance  amplitude  of  0.01. 

Even  at  M\  =  0.6  the  plots  of  uz  and  u>2/p  are  little  different  from  the  vorticity 
structure  for  the  incompressible  case,  found  by  Michalke  [1965a].  The  two  ‘elemen¬ 
tary  vortices’  in  the  eigensolution  will  subsequently  rotate  around  each  other  in  the 
non-linear  region  of  growth  and  merge  to  form  the  fundamental  vortex  in  the  mixing 
layer.  What  is  interesting  is  that  the  density  and  pressure  disturbance  fields  (figures 
2.22  c,d)  show  striking  similarities  with  the  fields  to  be  presented  in  Chapter  4  from 
the  non-linear  roll-ups  in  the  mixing  layer.  Low  density  and  pressure  perturbations 
are  found  in  the  vortex  core  and  high  density  and  pressure  are  found  in  the  region 
between  vortices  where  the  braid  will  eventually  form. 

Some  insight  into  mechanisms  can  be  obtained  by  examining  terms  in  the  com¬ 
pressible  vorticity  equations  for  the  linear  eigensolutions.  The  equations  for  uz  and 
for  w z/p  are  as  follows. 


du  dv 
Uz  dy  dx 


(2.55) 


D(jjz  _  / du  dv\  1  / dp  dp  dp  dp\ 
Dt  Uz\dx  +  dy)  +  p2\d  xdy  dydx) 

D{vz/p)  _  \_(  dp  dp  dpdp\ 

Dt  p3  \  dx  dy  dy  dx ) 


(2.55) 


(2.56) 
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The  terms  on  the  right  hand  side  are  plotted  on  figures  2.22 /  and  2.22 g.  Both  the 
baroclinic  and  dilatational  terms  are  negative  (dashed  contours)  in  precisely  those 
regions  where  the  elementary  vortices  develop,  and  act  to  reduce  the  vorticity  in  the 
region  where  vortex  roll-up  is  trying  to  occur,  Thus  both  these  terms  act  to  reduce 
the  growth  rate  of  the  two-dimensional  instability  as  Mach  number  is  increased. 
Recent  work  by  Lele  and  Shariff  (private  communication  [1989])  suggests  that  the 
advection  term  is  more  important  than  either  of  these  terms  and  is  the  main  reason 
for  the  stabilizing  effect  of  Mach  number. 

The  appearance  of  the  same  physical  processes  in  the  linear  eigensolutions  as 
are  observed  in  the  later  non-linear  development,  as  found  by  direct  numerical 
simulation,  may  help  to  explain  the  surprising  finding  that  the  linear  amplification 
rate  is  directly  proportional  to  mixing  layer  growth  rate. 


2.3  Linear  Instability  Model  for  the  Mixing  Layer 

Results  from  sections  2.2.1  and  2.2.5  suggest  that  linear  processes  may  be  im¬ 
portant,  even  in  the  fully  developed  mixing  layer.  It  is  therefore  worthwhile  to 
consider  a  simple  linear  instability  model  of  the  mixing  layer  flow.  Some  of  the  ar¬ 
guments  were  presented  by  Monkewitz  and  Huerre  [1982]  to  explain  their  successful 
prediction  of  shear  layer  growth  rate,  and  are  extended  here. 

2.3.1  Model 

Consider  first  the  incompressible  uniform  density  mixing  layer.  Experiments 
show  that  the  flow  is  dominated  by  the  primary  two-dimensional  instability,  and 
the  developed  structures  show  strong  coherence  in  the  spanwise  direction.  The  spa¬ 
tial  development  of  the  flow  can  be  described  by  the  successive  growth  of  linear 
instability  modes,  with  longer  and  longer  wavelengths.  The  phases  in  the  develop¬ 
ment  are  as  follows: 

(1)  Growth  of  the  fundamental,  most  unstable,  linear  instability  mode  of  the  initial 
profile.  This  mode,  frequency  u*F,  grows  until  the  layer  has  grown  by  a  factor 
of  two  in  width,  and  then  saturates  out,  neither  growing  nor  diminishing. 

(2)  Growth  of  the  first  subharmonic  wave,  frequency  w*F/2.  This  mode  grows 
and  eventually  results  in  pairing  of  two  of  the  fundamental  mode  structures. 
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After  the  pairing  the  layer  has  again  grown  by  a  factor  of  two  and  this  mode 
saturates. 

(3)  Growth  of  the  next  subharmonic  wave,  frequency  w^/4,  again  resulting  in 
pairing,  doubling  of  width,  and  saturation. 

The  essential  behavior  of  the  model  is  shown  on  figure  2.23.  The  exponential 
growth  in  a  measure  of  energy,  say  E,  of  each  mode  is  shown  schematically.  First 
the  fundamental  mode  (F)  grows  and  saturates,  and  we  associate  the  developed 
structure  with  the  neutral  instability  mode  (see  figure  2.24).  The  (non-dimensional) 
frequency  of  the  neutral  mode  o/jy  is  approximately  twice  the  frequency  of  the  most 
amplified  mode,  w,  where  w  is  non-dimensionalized  by  u>  =  u*6* /Uf.  Since  it  is 
the  same  wave  which  has  grown  and  saturated,  the  dimensional  frequency,  u>*  is  the 
same  for  the  two  cases,  and  the  layer  must  have  grown  by  a  factor  of  2.  Now  that 
the  layer  has  doubled  in  width  the  frequency  of  the  new  ‘most  amplified  wave’  is 
u>p/2,  the  first  subharmonic  wave  (Si).  This  grows  exponentially,  following  linear 
theory,  until  it  saturates,  resulting  in  pairing  of  two  of  the  neutral  modes  of  the 
original  instability.  The  resulting  structure  can  again  be  associated  with  a  neutral 
mode,  since  by  the  pairing  the  layer  has  grown  by  a  factor  of  2.  The  same  process 
of  successive  subharmonic  growth,  pairing,  and  saturation  at  the  neutral  mode, 
continues  ad  infinitum  in  the  streamwise  direction,  and  is  not  Reynolds  number 
dependent. 

The  essential  assumption  of  the  linear  instability  model  is  that  the  time  taken 
during  the  exponential  growth  of  an  unstable  wave  is  long  compared  to  the  time 
taken  for  the  ultimate  non-linear  process  of  pairing  and  saturation. 

^exponential  growth  ^non-linear  (2.58) 

However,  in  reality  both  linear  and  non-linear  processes  may  be  governed  by  similar 
physical  processes,  and  it  may  be  that  the  time  taken  in  the  two  stages  are  linked. 
So,  even  if  equation  (2.58)  is  not  satisfied  exactly,  the  linear  model  may  still  give 
good  predictions. 

Before  discussing  model  predictions,  we  note  that  in  reality  the  growth  of  sub¬ 
harmonic  waves  in  the  mixing  layer  is  dependent  upon  the  phase  of  the  growing 
wave  relative  to  the  large-scale  structure.  The  model  considers  only  linear  waves, 
for  which  phasing  is  not  important.  In  the  real  case  there  will  be  a  random  phasing 
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of  growing  waves,  which  has  the  effect  of  moving  the  pairing  locations  around  in 
space,  resulting  in  the  time  average  of  linear  growth  rather  than  as  a  series  of  steps, 
(see  also  Sandham  and  Reynolds  [1989]).  Modes  with  certain  phases  will  be  more 
amplified  by  the  presence  of  large-scale  structure,  and  modes  with  other  phases  will 
have  their  growth  rates  diminished  by  the  large-scale  structure.  The  model  relies 
on  these  effects  cancelling  out  in  the  long  time  average. 

2.3.2  Prediction  of  Growth  Rate 

The  first  implication  of  the  model  is  on  growth  rate.  Spatial  instability  waves 
grow  exponentially  like  e~aiX.  The  x  distance  for  the  n’th  most  amplified  wave  to 
grow  by  a  factor  of  eN  is  lV/|o^|max‘  hi  this  time  the  mixing  layer  has  doubled  in 
width,  so  AS  —  £n+1  —  Sn  =  Sn  or  A{S/Sn)  =  1,  so  the  growth  rate  is  given  by: 

A{S/6n)  1 

A{x/Sn )  lVy|Q!t'|max 

In  the  long  time  average  this  gives 

^  ~  |®i|max  (2.60) 

thus  providing  a  theoretical  basis  for  the  relation  found  in  section  2.2.1.  (equation 
(2.48)). 

It  should  be  noted  that  this  derivation  is  dependent  on  the  factor  N  being  a 
constant.  This  means  that  the  ratio  of  the  amplitude  of  a  wave  at  saturation  to  the 
amplitude  when  the  wave  was  the  most  amplified  wave  from  the  previous  saturated 
state  is  assumed  constant  for  all  x  locations. 

2.3.3  Prediction  of  Convective  Velocities 

A  second  corollary  of  the  model  is  that  the  large-scale  structures  found  in  the 
mixing  layer  are  linked  to  the  neutral  instability  modes.  The  model  picture,  figure 
2.23,  is  that  what  are  observed  in  the  mixing  layer  at  any  time  are  not  the  growing 
instability  waves,  but  the  neutral  mode  from  saturation  of  the  previous  instability. 
The  neutral  modes  are  steady  solutions  of  the  parallel  flow  equations  and,  since  they 
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are  neither  growing  nor  decaying,  would  persist  for  a  long  time  once  generated  in 
the  flow.  For  the  hyperbolic  tangent  velocity  profile  the  neutral  instability  modes 
are  a  subset  of  the  family  of  mixing  layer  vortices  found  by  Stuart  [1967],  The 
Stuart  vortices  were  compared  to  experimental  mixing  layer  structures  by  Browand 
and  Weidman  [1982],  The  comparison  was  qualitatively  good,  though  some  differ¬ 
ences  were  found.  Thus  it  may  be  that  the  neutral  instability  modes  are  a  way  to 
understand  features  of  the  organized  structures  in  free  shear  flows.  In  particular, 
it  is  postulated  here  that  the  convective  velocity  of  the  large-scale  structures  in  the 
flow  is  close  to  the  phase  speed  of  the  neutral  mode. 

It  was  found  that  the  phase  speeds  of  the  neutral  modes  were  dependent  on 
profile  shapes.  When  a  hyperbolic  tangent  was  used  for  the  mean  velocity  profile  it 
was  found  that  results  matched  the  convective  velocity  formula  (1.2)  almost  exactly. 
However,  when  the  solution  of  the  boundary-layer  equations  was  used  for  the  mean 
velocity  profile  it  was  found  that  the  phase  speed  of  the  neutral  mode  cpi  was 
different  to  the  Uc  formula  (1.2).  In  fact  it  was  always  biassed  towards  the  speed  of 
the  free-stream  with  the  highest  density,  as  shown  on  table  2.6  for  the  Brown  and 
Roshko  [1974]  velocity  and  density  ratios  (recall  1*2  =  I/P2),  at  =  0.1. 


Table  2.6  Phase  speeds  of  neutral  modes  (low  Mach  number) 


Mx 

U2 

t2 

CN 

Uc 

0.1 

0.378 

1.0 

0.712 

0.689 

0.1 

0.378 

0.143 

0.444 

0.549 

0.1 

0.378 

7.0 

0.943 

0.829 

0.1 

0.143 

1.0 

0.628 

0.572 

0.1 

0.143 

0.143 

0.243 

0.378 

0.1 

0.143 

7.0 

0.924 

0.765 

It  is  not  surprising  that  the  cjy  results  for  a  hyperbolic  tangent  mean  profile 
match  the  Uc  formula  (1.2).  The  hyperbolic  tangent  velocity  profile  has  an  anti¬ 
symmetry  about  the  centerline.  The  same  antisymmetry  is  in  the  structure  model 
from  which  the  convective  velocity  formula  (1.2)  is  derived.  This  effect  may  explain 
some  results  of  Lele  [1989],  who  performed  direct  two-dimensional  simulations  of 
spatially-developing  compressible  mixing  layers  beginning  with  a  hyperbolic  tangent 
mean  velocity  profile,  and  found  good  agreement  with  the  Uc  formula.  It  may  be 
that  this  observation  was  not  a  confirmation  of  the  Uc  formula,  but  rather  that  an 


38 


antisymmetry  was  built  into  the  whole  simulation  by  the  choice  of  inflow  velocity 
profile. 

Convective  speeds  for  the  compressible  mixing  layer  have  been  measured  ex¬ 
perimentally  by  Papamoschou  [1989].  Table  2.7  shows  a  comparison  between  the 
measured  convective  velocity,  the  prediction  of  the  Uc  formula  (1.2)  (matching  cjy 
results  from  a  hyperbolic  tangent  profile),  and  the  phase  speed  of  the  neutral  mode 
Cff  using  the  boundary-layer  mean  flow.  All  except  the  last  of  the  experimental 
results  show  phase  speed  skewed  towards  the  velocity  of  the  more  dense  stream,  as 
in  the  Cff  yielded  by  stability  theory  with  the  boundary-layer  mean  flow.  A  plot  of 
Mc\  versus  Mci  is  shown  on  figure  2.25.  This  kind  of  plot  was  used  by  Papamoschou 
to  demonstrate  how  different  the  actual  measured  convective  velocities  were  from 
the  Uc  prediction.  The  figure  shows  the  straight  line  Mc\  =  Mc 2  (which  is  the 
Uc  prediction),  Papamoschou’s  data  points,  and  the  phase  speeds  of  the  neutral 
modes,  cjy.  The  linear  model  seems  on  the  whole  to  do  a  better  job  of  predicting 
the  experimental  data  than  the  Uc  formula  (1.2). 


Table  2.7  Phase  speeds  of  neutral  modes  compared  with  Papamoschou  [1989] 


Mi 

U2 

t2 

cn 

Uc 

expt 

MEM 

mssm 

0.980 

0.978 

Erl 

0.906 

2.8 

BUI 

0.832 

0.857 

0.829 

1.7 

0.50 

0.109 

— 

0.512 

3.2 

0.13 

4.348 

0.903 

0.718 

0.878 

2.7 

0.13 

1.493 

0.771 

0.844 

2.6 

0.42 

0.182 

0.478 

0.593 

0.435 

3.1 

0.30 

0.400 

0.439 

0.571 

0.355 

3.2 

0.08 

1.205 

0.757 

0.561 

0.959 

3.0 

0.06 

0.535 

0.382 

0.457 

0.853 

These  results  led  Papamoschou  [1989]  to  propose  an  alternative  large-scale  struc¬ 
ture  with  shock  waves,  to  account  for  the  biassing  of  the  convective  velocity  towards 
one  or  other  of  the  free-streams.  The  contention  here  is  that  this  biassing  is  due  to 
the  experimental  mean  profiles  not  having  any  built-in  antisymmetry.  When  this 
effect  is  put  into  the  stability  analysis  (by  using  the  boundary-layer  mean  flow)  the 
biassing  is  captured  in  the  linear  theory,  without  resort  to  such  non-linear  effects 
as  shock  waves. 
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Other  experimental  measurements  of  convective  velocities  are  limited.  Brown 
and  Roshko  [1974J  measured  one  case,  with  U2  =  l/y/7  and  p2  =  7,  and  found 
a  structure  velocity  of  Uc  =  0.53,  compared  to  the  linear  estimate  of  0.444  and 
the  Uc  formula  estimate  of  0.55  (table  2.6).  For  this  case  u  was  found  that  the 
mean  profiles  of  velocity  and  density  plotted  by  Brown  and  Roshko  did  not  match 
well  with  the  solution  to  the  boundary-layer  equations.  We  note  that  the  solution 
for  mean  profiles  from  the  boundary-layer  equations  is  the  same  as  solving  the 
self-similar  mixing  layer  problem  for  the  mean  flow,  assuming  an  eddy-viscosity 
turbulence  model.  The  linear  theory  is  thus  limited  by  the  accuracy  of  the  mean 
flow  upon  which  it  is  performed. 


2.3.4  Prediction  of  Pairing  Locations 

Another  use  for  the  model  is  in  generalizing  the  findings  of  Bradshaw  [1966]  that 
approximately  1000  initial  momentum  thicknesses  are  required  in  the  streamwise 
direction  before  self-similarity  occurs,  and  of  Ho  and  Huerre  [1984]  who  modeled 
the  downstream  location  of  pairings  in  terms  of  the  initial  instability. 

Using  equation  2.486  above  and  the  linear  model  of  the  growth  process  (figure 
2.23),  we  can  write 


x  to  nth  pairing 


2n+1  -  1 
0.45 1  OC|  [max 


(2.61) 


where  x  is  normalized  as  x*/6*o  (6^  is  4  momentum  thicknesses  for  the  tanh  profile). 
The  n  =  0  event  is  the  roll-up  of  the  first  structure. 

As  an  example,  consider  the  mixing  layer  simulated  numerically  by  Lowery  and 
Reynolds  [1986].  This  was  an  incompressible  simulation,  with  equal  free-stream 
densities  and  a  velocity  ratio  U2  =  0.5.  From  table  2.2  we  have  that  |a^|max  is 
0.12850.  Using  equation  (2.61)  we  find  that  the  first  roll-up  is  predicted  to  occur 
at  x  =  17.3,  with  the  first  pairing  at  x  =  51.9  and  the  second  pairing  at  x  =  121.1, 
all  in  good  agreement  with  the  simulations.  The  Bradshaw  [1966]  criterion  for  self- 
similarity  would  occur  for  this  case  around  the  third  pairing  event  x  =  260.0,  or 
1040  in  intial  boundary-layer  momentum  thickness  units.  In  view  of  the  work  of 
Ho  et  al.  [1988],  who  found  phase  decorrelation  in  mixing  layers  to  be  associated 
with  the  pairing  events,  it  may  be  that  self-similarity  is  only  achieved  after  enough 


pairings  have  happened,  and  the  flow  has  forgotten  the  initial  conditions.  Three  or 
four  levels  of  pairing  would  seem  to  be  necessary  to  satisfy  Bradshaw’s  criterion. 

For  incompressible  mixing  layers  it  is  found  that  the  ratio  wneutrai/u;mo8t  amplified 
is  always  near  2.  However  for  compressible  mixing  layers,  using  the  boundary 
layer  velocity  profile,  this  was  not  found  to  be  always  true.  In  this  case  multiple 
interactions  ( e.g .  triplings)  may  be  more  common,  and  equation  (2.61)  can  be 
rewritten  as: 


x  to  nth  agglomeration  =  1  (2  62) 

0.45 1  a,- 1  max 

It  is  hoped  that  this  discussion  of  a  linear  instability  model  for  the  mixing  layer 
prompts  more  experimental  measurements  of  growth  rate,  mean  profiles,  and  the 
details  of  the  large  scale  structure  in  the  flow,  including  convective  velocities,  average 
pairing  locations,  number  of  structures  involved  in  agglomerations  etc.  In  this  way 
the  limitations  of  the  model  can  be  explored  and  improvements  made. 

2.4  Chapter  Summary 

In  this  chapter  it  has  been  shown  that  the  linear  theory  can  be  a  very  useful  tool 
for  investigating  the  compressible  mixing  layer.  In  particular  : 

(1)  The  maximum  amplification  rate  found  from  spatial  theory,  using  the  solution 
of  the  compressible  boundary- layer  equations  as  the  mean  flow,  appears  to 
be  directly  proportional  to  the  growth  rate  of  the  developed  mixing  layer 
(equation  (2.48)).  This  was  demonstrated  for  the  effect  of  density  and  velocity 
ratios  on  mixing  layer  growth  rate  at  low  Mach  number. 

(2)  Oblique  waves  are  dominant  in  the  mixing  layer  above  a  convective  Mach 
number  of  0.6.  The  obliquity  of  the  most  amplified  wave  selected  by  the  linear 
theory  is  such  that  the  convective  Mach  number  perpendicular  to  the  wave 
crest  is  approximately  0.6  (equation  (2.53)). 

(3)  The  experimental  observation  of  Papamoschou  [1989],  that  actual  convective 
velocities  did  not  agree  with  the  Uc  theory  (equation  (1.2)),  was  confirmed 
in  the  linear  regime  by  calculating  the  phase  speeds  of  the  neutral  modes  of 
instability. 
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CHAPTER  3 


Numerical  Formulation  for  Direct  Simulations 


In  a  direct  simulation  the  time-dependent  Navier-Stokes  equations  are  solved  in 
full,  with  no  turbulence  model.  This  chapter  describes  the  numerical  methods  used 
to  simulate  the  compressible  time-developing  mixing  layer,  which  is  assumed  to  be 
periodic  in  the  streamwise  and  spanwise  directions  (x  and  z).  The  code  was  written 
in  the  VECTORAL  language  (Wray  [1988]),  and  implemented  on  a  Cray  X-MP  4/8 
at  NASA-Ames. 


3.1  Governing  Equations 


The  conservation  laws  for  mass,  momentum  and  energy  can  be  written  in  the 
following  form  (Anderson,  Tannehill  and  Pletcher  [1984]),  using  the  conservative 
variables  (p* ,  p*  u* ,  E^) ,  cartesian  tensor  notation,  and  the  superscript  *  for  a  di¬ 
mensional  quantity: 


dp*  d(p*u*) 
dt *  dx* 


=  0 


(3.1) 


d{p*u*)  d(p*u*u*  +p*6X])  dr*3 

dt *  +  dx*  ~  dx* 

3  3 


(3.2) 


dE’r  a|(gf  +  )■>,•]  = 

dt*  dx*  dx*  dx* 

l  «  » 


(3.3) 


where  p*  is  the  density,  u?  are  the  velocity  components,  p*  is  the  pressure,  r*.  is  the 
shear  stress  tensor  and  q*  is  the  heat  flux  vector.  Ej  is  the  total  energy,  defined 
by: 

E't  =  />*(«*  +  (3-4) 

where  e*  is  the  internal  energy  per  unit  mass.  The  perfect  gas  law  is 

p*  =  p*R*T*  (3.5) 


43 


where  R*  is  the  gas  constant  and  T*  the  temperature.  We  assume  a  Newtonian 
fluid  and  Fourier  heat  conduction,  so  the  constitutive  equations  for  t*-  and  q*  are: 


_*  =u*(dUi  |  2dUkS  \ 

'i  M  \dx*  +  8x*  3dx*StJ) 


(3.6) 


(3.7) 


where  n*  is  the  viscosity  and  k*  is  the  thermal  conductivity. 
Non-dimensionalization  of  these  equations  is  obtained  by: 


u* 

t 

•1  .  -  * 

_£l 

p * 

t-TL 

(3.8a) 

* “  v; 

P  P\ 

P  n*TT*2 

Ti 

(*• 

M  =  — 

e * 

e  =  —T 

Uf 

t*u: 

1  - «. 

H 

«*. 

II 

(3.8ft) 

where  the  subscript  1  represents  the  upper  (y  >  0)  free-stream  value.  The  reference 
length  scale  is  6*o,  the  vorticity  thickness  of  the  initial  velocity  profile: 


vi  -  u2 

\du*/dy*\  max 


(3.9) 


The  non-dimensional  equations  for  continuity,  momentum  and  energy  are: 


dp  _  d{puj) 

dt  dxi 


d{puj)  =  d(puiUj  +  p6jj) 
dt  dxj 


dEf  _  9[(Et  +  p) tzt- ]  dq 

dt  dx%  dii 


(3.10) 

drjj 

dxj 

(3.11) 

d(UjTij) 

dii 

(3.12) 
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with  constitutive  relations 


=  _M_  (dui  i  duJ  2duk , 

i2e\dx,  dxt  3  3x*. 


(3.13) 


(^f  —  l)MjPr/2e  dz» 


(3.14) 


The  Reynolds  number  of  the  flow  is  defined  by  Re  =  p\U{t>w0l A*i  >  and  the  Prandtl 
number  by  Pr  —  c*p* /k*.  The  Prandtl  number  is  assumed  constant  with  value  1. 
The  viscosity  is  assumed  to  follow  a  power  law,  so  non-dimensionally 

H  =  T0-67  (3.15) 

The  constant  0.67  is  the  value  for  nitrogen  (White  [1974]). 

Constant  specific  heats  are  assumed,  and  if  we  set  e*  =  0  at  T*  —  0,  we  can  write 
e*  =  c„T*.  The  non-dimensional  form  of  the  perfect  gas  law  can  then  be  written  as 
either: 

P  =  p[~1-  l)e  (3.16) 


IP  =  PC2  (3.17) 

where  c  is  the  sound  speed  (normalized  as  c* /U^)  and  we  have  used  the  result  for 
the  sound  speed  of  a  perfect  gas  with  constant  specific  heats: 


c*  =7i?*r 


(3.18) 


In  all  the  simulations  a  passive  scalar  field  is  computed.  The  dimensional  equation 
for  the  scalar  /  (see  e.g.  Kays  and  Crawford  [1980])  is 

a(p'i)  .  9(p 7<)  a 


at 

where  D*  is  the  diffusion  coefficient  for  the  scalar.  One  can  think  of  p*  f  as  the 
concentration  per  unit  volume  of  a  trace  species.  We  define  a  Schmidt  number  by 
Sc  =  p*  jp*D\  and  using  the  non-dimensionalization  method  of  equation  (3.8)  we 
have: 

d{pf)  djpfuj)  =  1  3  (  df 

dt  dx .  ReScdx 


3.2  Time  Advance 


The  time  advance  method  is  fully  explicit.  The  variables  [p,pu,pv,pw,Ej)  are 
advanced  using  a  three-step  compact-storage  third  order  Runge-Kutta  scheme  of 
the  family  derived  by  Wray  [1986].  Two  storage  locations  are  employed  for  each 
time-dependent  variable  and  at  each  substep  these  locations,  say  Qi  and  Q2,  are 
updated  simultaneously  (using  a  feature  of  VECTORAL)  as  follows: 

Qf w  =  aiAfQf1  +  Qfd  Q2ew  =  a2AfQ?ld  +  Q|ld  (3.21) 

The  constants  (ai,a2)  are  chosen  to  be  (2/3, 1/4)  for  substep  1,  (5/12,3/20)  for 
substep  2  and  (3/5, 3/5)  for  substep  3.  At  the  beginning  of  each  full  time  step 
Ql  and  Q2  are  equal.  The  data  in  Qi  is  used  to  compute  the  right  hand  side  of 
equations  (3.10)  through  (3.12).  The  computed  right  hand  side  is  stored  m  Qi 
(overwriting  the  old  Q\).  Equation  (3.21)  is  then  used  to  update  Q\  and  Q2. 

An  estimate  of  the  time  step  limitation  can  be  made  by  considering  a  model 
convection-diffusion  equation  (Blaisdell  [1988],  private  communication): 

d<f>  d<f>  d<j>  d<f>  d*<j> 

—  +  u— — h  v— — h  xv—  =  a— — - — 
dt  dx  dy  dz  dx^dx^ 

Assuming  periodic  boundary  condition  in  all  directions  we  can  Fourier  transform 
and  rearrange  this  equation  to  give: 


~l7) 


+  a47r2 


+ 


a 


(3.23) 


This  equation  can  be  compared  with  a  simple  ordinary  differential  equation: 


for  which  the  time  step  limitation  for  stability  is 


{CF  L)  majc 

|&|max 


(3.24) 


(3.25) 
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where  CFL  is  the  Courant-Friedrichs-Lewy  number.  The  maximum  CFL  number 
for  stability  is  fixed  by  the  time-advance  method.  For  linear  equations  the  limit  is 
y/Z  for  the  third  order  Runge-Kutta  method  described  above. 

The  magnitude  of  a  complex  number  a+ib  is  less  than  |a|-|-|6|,  so  we  can  substitute 
absolute  values  for  the  terms  on  the  right  hand  side  of  equation  (3.23).  We  choose 
the  worst  case  for  the  terms  in  brackets  in  (3.23),  ».e.  kx  =  Nx/ 2,  u  =  c  +  |u|,  etc. 
The  value  for  a  in  (3.22)  is  obtained  from  (3.14)  as  n/fa  —  1  )M^RePr,  which  is  a 
more  stringent  limitation  than  the  momentum  equation  for  <1.6  (for  Pr  =  1 
and  7  =  1.4).  Letting  Ax  =  Lx/Nz  and  similarly  for  Ay  and  A z  we  get: 


where 


CFL 

Dc  +  Dfi 


Dc  =  7TC 


/  1  1  1  \  /  |u|  |v|  ju; 

vAx  +  Ay  +  A  z)  +  "Ux  +  Ay  +  Az 

*2t*/p  (  1  JL  ,  1  \ 

(7  -  l)M%RePr  ''Ax2  Ay2  A z2) 


(3.26a) 


(3.266) 


(3.26c) 


The  worst-case  cell  is  used  to  fix  the  time  step. 

Simulations  with  various  time  steps  were  made  for  a  time-developing  mixing  layer 
problem  at  Afj  =  0.4  and  Re  =  400,  initialized  with  the  fundamental  mode  from 
linear  stability  analysis.  The  simulations  were  run  up  to  the  non-linear  saturation  of 
the  instability.  It  was  found  that  there  was  no  numerical  instability  or  degradation 
of  accuracy  until  the  CFL  number  was  raised  above  4,  indicating  that  the  above 
criteria  is  very  conservative.  In  all  the  simulations  a  CFL  number  of  2  was  used. 


3.3  Evaluation  of  Derivatives 
3.3.1  Periodic  Directions 

In  the  periodic  directions  x  and  z  the  Fourier  method  is  used  to  calculate  deriva¬ 
tives.  The  quantity  to  be  differentiated  is  transformed  using  the  Fast  Fourier  Trans¬ 
form,  multiplied  by  »'fc,  and  inverse  transformed.  No  attempt  was  made  to  dealias 
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the  non-linear  products,  since  no  complete  method  is  known  for  the  compressible 
equations  which  have  many  terms  involving  divisions. 

When  an  even-length  Fourier  transform  is  used  there  is  one  wavenumber  kx  = 
Nx/ 2  which  has  no  complex  part.  This  wavenumber  is  commonly  referred  to  as 
the  ‘oddball’  wavenumber.  The  Fourier  component  at  the  oddball  wavenumber 
is  always  zeroed  out  when  taking  derivatives.  In  an  implicit  scheme  this  Fourier 
component  is  forced  to  remain  zero.  However  in  am  explicit  scheme,  such  as  that 
used  here,  roundoff  error  will  accumulate  at  this  wavenumber  as  one  goes  back 
and  forth  between  Fourier  space  and  physical  space.  In  the  present  code  this  is 
compounded  by  conversions  back  and  forth  between  64  bit  arithmetic  representation 
and  32  bit  storage  representation  (Blaisdell  [1988],  private  communication).  It  is 
thus  necessary  to  explicitly  zero  out  the  oddball  component  of  the  flow  during  any 
computation.  Blaisdell’s  procedure  for  doing  this  (for  one  direction)  is  to  define  a 
filtered  flowfield  <f>j ,  subtracting  the  oddball  as  follows: 

</,'  =  (3.27) 

From  the  definition  of  a  transform  pair,  we  have 

N ,/2  u 

*j  =  £  V  (3.28) 

kx=  —  -jj*-  + 1 

?(*«)  =  Y.  (3.29) 

J  =  1 

Substituting  for  <7^/2  'n^°  equation  3.27  gives 

<3-3o> 

1  i 
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This  method  of  removing  the  oddball  in  real  space  vectorizes  efficiently  and  was 
found  to  be  substantially  cheaper  than  the  alternative  method  of  Fourier  trans¬ 
forming,  zeroing  the  oddball,  and  inverse  transforming.  The  method  was  applied 
at  each  substep  of  the  integration. 

In  the  three-dimensional  code  the  oddball  needs  to  be  removed  in  each  spectral 
direction,  z  and  z.  Following  the  same  procedure  as  above,  the  filter  function  is 
derived  as: 


(3.31) 


,f  i  (-1)fc  v'/  *\k'i  (_1)J  ,w'i 

+j,k  =  tj*  ~  Is  ~  “atT  Is  (-1)  W 


k'= 1 


j'= 1 


(— 1)J (~ l)fc  Y' ,  .xi*/ 

+  ^  ^ (-1)  (_1) 


x^z 


j'=ik'=i 


3.3.2  Normal  Direction 

In  the  normal  direction  ( y )  there  are  several  alternative  differentiation  methods 
which  can  be  used.  Firstly,  one  can  subtract  the  mean  flow  from  the  quantity 
to  be  differentiated,  and  use  regular  Fourier  differentiation.  This  would  only  be 
appropriate  for  a  periodic  array  of  mixing  layers,  which  is  not  the  case  we  wish 
to  consider  here.  Secondly,  one  can  use  a  mapping  ( e.g .  Cain  ct  al.  [1984]),  to 
transform  an  infinite  physical  domain  onto  a  finite  computational  one.  This  method 
was  implemented,  but  was  found  to  produce  point-to-point  oscillations  in  the  free- 
stream.  This  effect  is  probably  due  to  a  lack  of  resolution  of  sound  waves  in  the 
free-stream,  where  the  grid  spacing  becomes  very  large.  The  method  finally  adopted 
was  to  use  high-order  finite  difference  methods  on  a  grid  of  finite  extent  in  the 
normal  direction,  and  then  to  apply  characteristic  boundary  conditions  to  simulate 
an  infinite  flow  domain. 

A  family  of  high  order  modified  Pade  schemes  has  been  derived  by  Lele  [1989]. 
These  schemes  can  be  written  in  the  following  form  for  the  first  derivative: 

ii  ,  .ii  ,  ii  .  .  >j'+2  - 

+  a4>i  +  +i- 1  -  b - 5a7 - +  c—  7ZZ  (3-32a) 
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Solution  for  <f>'^  is  obtained  by  solving  the  tridiagonal  system  of  equations.  A  family 
of  fourth  order  schemes  are  obtained  if  we  choose: 


2  +  4  a 
3 


c  = 


4  —  o 
3 


(3.326) 


For  o  =  4,  the  conventional  Pad£  scheme  is  recovered,  whilst  with  a  =  3  we  have  a 
sixth  order  scheme,  used  in  all  the  current  work.  A  similar  scheme  can  be  derived 
for  the  second  derivative: 


r  ,+af+r  =  z ty+ii-i 

Vj-i  -ru<Pj  TVj-i  o  Aj/2  -rc  4Ajj2 


4A  y2 


(3.33a) 


,  4a  —  4  10  —  a 

b  = -  c  - - 

3  3 

This  time  a  =  10  gives  the  usual  Pade  scheme,  and  a  =  5.5  is  a  sixth  order  scheme, 
used  in  this  work. 

The  viscous  terms  in  the  governing  equations  require  evaluation  of  successive 
derivatives,  for  example  the  term 


d  (  du\  .  . 

(3-34) 

When  a  spectral  method  is  used  there  is  no  loss  of  accuracy  if  these  are  computed 
by  two  applications  of  a  first  derivative.  However,  with  finite  difference  methods  we 
find  that  two  applications  of  a  first  derivative  gives  a  much  worse  representation  of 
the  high  wavenumbers  than  a  single  second  derivative  computation.  This  is  because 
the  modified  wavenumber  (Lele  [1989])  goes  to  zero  for  the  first  derivative  at  high 
wavenumbers.  The  remedy  for  this  is  to  expand  all  terms  in  the  y  direction  like 
(3.34)  into  two  terms  (non-conservative  formulation): 


d2u  d/idu 
^  dy2  +  dy  dy 


(3.35) 


A  mapped  grid  is  used  in  the  normal  direction  to  resolve  the  shear  layer  more 
efficiently.  The  mapping  function  is  a  sinh  function,  which  concentrates  points 
around  y  =  0  (Anderson,  Tannehill  and  Pletcher  [1984]): 


Ly  sinh(6yT7) 
2  sinh  by 


(3.36) 
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where  Ly  is  the  box  length  in  the  normal  direction.  The  mapped  coordinate  is  rj, 
and  by  is  the  stretching  parameter  in  the  y  direction,  chosen  to  be  1.7  in  the  current 
work.  If  we  define  the  metrics 


h,  _  dy  _  Ly  by  cosh(feyq) 
dr)  2  sinh  by 


(3.37) 


„  =  =  Ly  fey  sinh(6y>?) 

dr)2  2  sinh6y 

then  the  first  and  second  derivatives  of  a  function  4>  are  evaluated  as: 


(3.38) 


d<f>  _  1  d<f> 
dy  h!  dr) 


(3.39) 


d2<f>  _  1  d2y  hH  d<f> 
dy2  h 12  dr)2  hf2  dy 


(3.40) 


3.4  Free-Stream  Boundary  Condition 

In  these  simulations  we  are  considering  the  problem  of  an  unbounded  compress¬ 
ible  mixing  layer.  The  infinite  extent  could  be  obtained  by  using  a  mapping  function, 
but  this  would  lead  to  poor  resolution  of  the  flow  far  away  from  the  mixing  layer. 
In  particular,  sound  waves  would  propagate  into  regions  of  the  computation  where 
they  would  be  poorly  resolved,  and  might  be  reflected  back  and  contaminate  the 
main  flow.  Thus,  we  require  boundary  conditions  which  simulate  an  infinite  domain, 
even  though  the  computational  domain  is  finite.  Such  schemes  were  investigated  by 
Thompson  [1987],  and  his  method  is  followed  here.  The  basic  idea  is  to  consider  the 
characteristic  form  of  the  Euler  equations  at  the  boundary.  Outgoing  characteris¬ 
tics  use  information  within  the  computational  domain,  and  can  be  computed  with 
no  problem.  Incoming  characteristics  are  handled  by  setting  the  time  derivative  of 
their  amplitude  equal  to  zero,  thus  ensuring  that  no  waves  enter  the  domain  during 
the  simulation,  giving  the  boundary  conditions  a  non-reflecting  character. 


We  begin  by  writing  the  Euler  equations  in  terms  of  the  conservative  variables 
Q  =  [p,Pu,pv,pw,ET,pf)T 


dQ  ,  dG 

ST  +  -Ty  =  (r-h-s) 


(3.41) 


G  =  ( pu,puv,pv 2  +  ptpvwt  (Ey  +  p)v,pvf)T  (3.42) 

We  are  concerned  here  only  with  the  Euler  derivatives  in  the  y  direction.  Derivatives 
in  x  and  z,  including  viscous  terms,  are  evaluated  at  the  boundary  using  informa¬ 
tion  from  the  previous  substep,  and  are  included  in  the  right  hand  side  (r.h.s).  The 
flux  Jacobian  of  G  is  more  easily  derived  if  we  work  with  the  non-conservative 
flow  variables  U  =  (p,u,t/,u;,d,  f)T,  where  d  =  p/»“7.  Setting  B  —  dG/dQ 
(».e.  J3(row  t,col  j)  =  dG(i)/dQ(j))  and  R  =  dQ/dU,  we  have 


and 


au  nau 

dt  +B  dy 


R  1(r.h.s) 


(3.43) 


dG  au 

dy  ~RB  dy 


(3.44) 


Now  B  can  be  diagonalized,  B  =  T-1A!T,  where  the  elements  of  the  diagonal  matrix 
A  are  (v  —  c,v,v,v,v  +  c,v).  Pressure  is  more  easily  found  computationally  than 
the  quantity  d,  so  we  use  a  new  non-conservative  vector  V  =  (p,u,v,w,p,  f)T  and 
define  S  =  dU/dV.  Equation  (3.44)  can  now  be  written 


dG 

dy 


(3.45) 


This  is  the  relation  that  is  used  to  calculate  dG/dy  in  equation  (3.41)  at  the  bound¬ 
ary. 

The  quantity  in  brackets  in  equation  (3.45)  is  a  vector.  The  sign  of  each  of  the 
eigenvalues  in  A  is  used  to  determine  the  course  of  action  for  each  element  in  the 
vector.  If  the  characteristic  velocity  is  directed  out  of  the  computational  domain 
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(positive  at  the  top  boundary  or  negative  at  the  bottom  boundary),  then  the  quan¬ 
tity  is  calculated  as  is.  On  the  other  hand,  if  the  characteristic  is  directed  inwards 
then  the  vector  element  is  zeroed  out.  This  gives  the  non-reflecting  character  of 
the  boundary  condition.  The  final  step  is  to  premultiply  the  vector  by  the  matrices 
T'1  and  R.  The  various  vectors  and  matrices  required  in  the  computation  axe: 


A  TS—  = 
dy 


.,dw 

<«+«)if5 +<*$*) 

vU 

dy 


(3.46) 
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(3.48) 


c2  u2  -(-  V2  +  w 2 

where  a  = - 1 - - - 

7-1  2 


The  boundary  condition  for  the  two-dimensional  code  is  a  simple  extract  from 
the  above,  removing  the  fourth  row  and  column  from  each  vector  and  matrix. 
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3.5  Initial  Conditions 


To  start  the  calculation  we  need  to  specify  the  values  of  all  computational  vari¬ 
ables.  The  specifications  are  described  in  subsequent  chapters. 


3.6  Validation 

The  objective  is  to  subject  the  code  to  as  many  checks  as  possible.  We  rely  on 
three  methods:  (a)  validation  in  the  incompressible  limit  using  Taylor-Green  flow, 
(b)  checking  the  Euler  terms  in  the  compressible  equations  by  comparison  of  linear 
growth  rates  with  inviscid  linear  stability  theory,  and  (c)  checks  of  resolution  by 
monitoring  energy  spectra  during  simulations,  to  ensure  that  enough  modes  are 
being  used  to  fully  resolve  the  flow.  The  only  terms  that  cannot  be  verified  using 
these  methods  are  the  viscous  terms  involving  du^/dx^  (equation  (3.13)).  Here  we 
must  rely  on  thorough  sight-checking  of  the  code. 

The  Taylor-Green  flow  (Taylor  [1923])  consists  of  a  decaying  array  of  vortices  in 
the  x  —  z  plane,  specified  non-dimensionally  by: 

u  =  -\f~A  cos(2ickx)  sm(2irkz)e~8r2k2t^Re  (3.49) 

w  =  \/Asin(27rfci)  cos(27rfcz)e~8,r  k  t/Re  (3.50) 

P  ~  Pref  —  —  [cos(47rfci)  ~t-  cos(47rfc2)]e-8,r  k  t!Re  (3.51) 

4 

We  choose  k  =  1,  Re  =  1,  a  box  size  Lx  =  Lz  =  1,  and  uniform  conditions  in  y  with 
Ly  =  2.  The  reference  pressure  is  pref  =  1  / (7M2).  The  reference  Mach  number 
Mi  was  set  to  0.2  and  the  constant  A  was  ta;.en  a s  0.0016,  so  that  we  begin  with  a 
peak  local  Mach  number  of  approximately  0.01,  i.e.  near  the  incompressible  limit. 
The  Taylor-Green  problem  needs  only  a  few  Fourier  modes.  The  test  case  was  run 
for  the  three-dimensional  code  with  Nx  x  Ny  x  Nz  =  8  x  11  x  8  corresponding  to 
4  complex  Fourier  modes  in  each  of  the  x  and  z  directions.  The  solution  for  u  at 
x  =  0.5  and  z  —  0.25  is  shown  on  figure  3.1  as  a  function  of  time.  The  error  at 
time  t  =  l/(87r2),  when  the  solution  should  have  decayed  by  a  factor  e,  is  less  than 
0.03%.  This  remaining  error  is  attributed  to  the  slight  difference  in  Mach  number. 
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The  growth  of  linear  disturbances  in  the  three-dimensional  code  was  checked 
against  the  inviscid  linear  stability  theory.  The  code  was  initialized  with  the  eigen¬ 
functions,  and  run  at  high  Reynolds  number  through  the  linear  regime.  The  linear 
growth  rate  was  extracted  from  the  Fourier  coefficients  from  the  relation 

d 

u(ks)  =  —  ln{Real[^(fcz)]}  (3.52) 

at 

For  the  fundamental  mode  we  have  kx  =  1.  Plots  of  ln(Real[<£(l)])  are  shown  on 
figure  3.2  for  the  functions  <f>  —  (u,v,w, p,T,p).  The  simulation  was  for  a  45° 
oblique  wave  growing  in  a  mixing  layer  at  Mi  =  0.8,  Re  =  106.  All  functions 
show  the  correct  linear  growth.  In  addition,  the  eigenfunctions  at  the  end  of  the 
simulation  were  renormalized  and  compared  with  the  initial  eigenfunctions.  An 
example  is  shown  on  figure  3.3  for  temperature.  All  the  eigenfunctions  were  found 
to  collapse  on  top  of  each  other.  The  worst  error  in  the  growth  rate  compared  to  the 
inviscid  theory  was  approximately  2%,  attributed  to  the  finite  Reynolds  number, 
and  changes  in  the  mean  flow  (setting  up  a  small  tJ  velocity  profile)  during  the 
simulation.  The  latter  effect  could  be  removed  by  running  for  just  one  time  step, 
when  the  mean  flow  had  no  time  to  change;  in  that  case  the  error  dropped  to  0.01%. 

The  effect  of  the  positioning  of  the  boundaries  in  the  y  direction  was  investigated 
by  running  three  simulations  at  M\  =  0.4,  Re  =  400.  The  length  of  the  domain 
in  the  streamwise  direction  was  Lx  =  8.0554  and  the  length  in  the  y  direction  was 
varied,  with  Ly  set  to  6,  10  or  20.  Ly  =  10  means  that  the  domain  extends  to 
y  =  ±5  initial  vorticity  thicknesses  on  either  side.  The  number  of  grid  points  in  y 
was  chosen  to  be  61,  81  and  99  respectively.  The  growth  of  the  vorticity  thickness 
is  plotted  on  figure  3.4  for  each  case,  and  the  developed  structure  at  the  same  time 
in  the  simulation  is  shown  on  figure  3.5  for  Ly  =  10  and  figure  3.6  for  Ly  =  6.  The 
physics  of  such  simulations  are  discussed  in  more  detail  in  the  following  chapter. 
Here  we  note  that  with  Ly  >  10  we  get  a  converged  solution  for  the  growth  history  of 
the  flow.  However,  even  the  case  with  Ly  =  6  captures  most  of  the  interesting  roll¬ 
up  of  large-scale  structure.  The  detailed  feature  which  arises  when  the  boundaries 
are  too  close  is  the  development  of  low-level  vorticity  at  the  boundary,  shown  in 
figure  3.66.  All  calculations  reported  herein  were  run  with  Ly  >  10. 
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One  of  the  most  important  checks  on  the  code  is  that  of  resolution.  In  general  if 
the  same  result  is  obtained  with  two  grids,  one  with  twice  as  many  points  in  each 
direction  as  the  other,  then  we  have  some  confidence  in  the  result.  With  a  spectral 
method  there  is  an  inbuilt  check  on  resolution,  namely  the  fall-off  in  the  energy 
spectrum  at  high  wavenumbers.  A  large  drop-off  in  energy  at  high  wavenumbers 
means  that  the  flow  is  fully  resolved,  and  adding  modes  would  not  improve  the 
solution.  Any  lack  of  resolution  will  show  up  as  an  upturning  in  the  spectrum  at 
high  kx •  During  the  simulations  the  energy  spectrum  is  monitored.  The  energy  is 
defined  by: 

rLy/2 

E{kx)  —  I  Ui(kx)uAkx)dy  (3.53) 

J-Ly/2 

where  t  indicates  a  complex  conjugate.  A  typical  well-resolved  spectrum  for  the 
Ly  =  10  case  above  is  plotted  on  figure  3.7,  showing  9  orders  of  magnitude  roll-off 
in  energy.  Monitoring  of  the  spectrum  is  also  used  to  improve  the  efficiency  of  the 
code.  When  linear  modes  are  used  for  initial  conditions  the  code  needs  few  modes 
to  resolve  the  flow.  As  the  harmonics  develop,  and  the  flow  becomes  non-linear, 
we  need  more  and  more  modes  to  resolve  the  flow,  until  the  flow  is  fully-developed 
and  the  small  diffusion  scales  get  no  smaller.  The  simulations  were  usually  started 
with  few  modes,  e.g.  16  x  99  x  16  and  then  as  the  spectrum  filled  out  the  spectral 
directions  x  and  z  were  extended,  ending  up,  for  example,  at  96  x  99  x  96.  In  the  y 
direction  we  need  the  full  number  of  points  right  from  the  start  in  order  to  resolve 
the  mean  flow.  Factors  of  2  were  typically  saved  in  run-time  using  this  method. 

As  a  last  point  on  verification  it  is  noted  that  the  two-dimensicral  simulations 
from  the  current  work  have  been  compared  with  results  from  an  earlier  code  (Sand- 
ham  and  Yee  [1989]),  which  used  TVD  and  MacCormack  methods,  and  with  the 
simulations  of  Lele  [1989].  All  such  comparisons  show  qualitatively  similar  behavior. 
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CHAPTER  4 

Two-Dimensional  Simulations 


In  this  chapter  results  are  presented  from  direct  simulations  of  the  compressible 
Navier-Stokes  equations  in  two  dimensions.  The  simulations  show  the  effect  of  com¬ 
pressibility  on  the  development  of  the  primary  two-dimensional  instability,  which 
has  been  observed  to  be  the  dominant  mode  in  incompressible  experiments. 


4.1  Initial  Conditions  and  Parameters 

The  simulations  have  been  performed  for  the  temporal  development  of  the  two- 
dimensional  time-developing  mixing  layer.  Mean  profiles  of  velocity  and  tempera¬ 
ture  need  to  be  specified  at  the  initial  time.  The  mean  velocity  in  the  streamwise 
direction  is  given  by  an  error-function: 

u  =  er^yv/^)  (4-1) 

The  mean  temperature  is  obtained  from  the  mean  velocity  profile  via  the  Crocco- 
Busemann  relation  (equation  (2.13)),  which  assumes  parallel  flow  and  unity  Prandtl 
number.  Pressure  is  assumed  uniform,  and  then  density  is  obtained  simply  as  the 
inverse  of  temperature,  since  both  density  and  temperature  are  normalized  by  the 
value  on  the  upper  side  of  the  mixing  layer  (equation  (3.8)). 

Perturbations  are  added  to  the  mean  profiles  in  the  form  of  eigenfunctions  of 
unstable  modes  from  the  inviscid  linear  stability  analysis  (Chapter  2).  For  example: 

u  =  u  +  A  Real{u(y)e,ai}  (4.2) 

and  similarly  for  the  v  velocity,  density  and  temperature.  The  amplitude  A  of 
the  disturbances  was  usually  specified  to  be  0.05.  The  length  of  the  box  in  the 
streamwise  direction  ( Lx )  was  chosen  to  match  the  most  unstable  wavelength  of 
inviscid  linear  stability  theory  i.e.  Lx  =  27r/|u^|max- 

Non-dimensional  parameters  for  the  flow  also  have  to  be  specified.  For  most  of 
the  simulations  in  this  chapter  an  initial  Reynolds  number  Re  =  400  was  chosen 
based  on  the  following  considerations.  The  Reynolds  number  has  to  be  chosen 
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small  enough  so  that  the  flow  can  be  fully  resolved,  yet  high  enough  to  capture 
the  inviscid  nature  of  the  instability.  The  effect  of  Reynolds  number  on  the  linear 
growth  rate  of  the  most  amplified  inviscid  eigenfunction  was  investigated  by  running 
the  viscous  code,  initialized  with  an  inviscid  eigenfunction,  for  one  time  step.  The 
linear  growth  rate  can  be  found  directly  from  the  change  in  the  Fourier  coefficient 
for  the  relevant  mode  using  equation  (3.52).  Results  are  shown  in  figure  4.1.  The 
asymptotic  growth  rate  at  high  Reynolds  numbers  matches  the  inviscid  theory  to 
within  0.1  percent.  The  drop  in  amplification  as  viscosity  is  increased  can  be  clearly 
observed.  It  should  be  noted  that  only  the  growth  rate  of  the  inviscid  eigenfunction 
is  plotted.  There  will  be  a  viscous  eigenfunction,  with  a  longer  wavelength,  that  is 
more  amplified.  Thus  the  curves  at  low  Reynolds  numbers  for  the  most  amplified 
viscous  wave  will  not  drop  so  rapidly  as  these.  However,  these  plots  do  give  a  good 
idea  of  the  range  of  Reynolds  numbers  where  we  expect  to  capture  the  basic  inviscid 
processes,  as  might  be  found  in  high  Reynolds  number  experiments.  In  choosing 
Reynolds  numbers  for  the  simulations,  we  would  like  at  least  to  be  in  the  region 
where  the  inviscid  eigenfunction  is  linearly  amplified,  and  not  damped.  It  can  be 
seen  from  the  plot  that  as  the  Mach  number  is  increased  we  have  to  move  to  higher 
and  higher  Reynolds  numbers  to  get  near  the  asymptotic  high-Re  region. 

The  effect  of  Reynolds  number  on  the  non-linear  development  of  a  single  struc¬ 
ture  at  Mi  =  0.4  was  computed.  The  growth  of  vorticity  thickness  is  plotted  on 
figure  4.2,  and  the  scalar  and  vorticity  fields  for  the  Re  =  100  and  Re  =  800  cases 
are  shown  on  figures  4.3  and  4.4.  At  even  lower  Reynolds  numbers,  with  the  same 
low  amplitude  forcing,  the  layer  is  laminar  and  grows  with  thickness  proportional 
to  the  square  root  of  time.  As  the  Reynolds  number  is  increased,  it  can  be  seen  that 
the  width  of  the  strained  diffusion  layer  in  the  stagnation  region  (x  =  3Lz/4,  y  =  0) 
is  reduced.  Broadwell  and  Breidenthal  [1982]  predict  that  the  width  of  this  layer 
varies  in  inverse  proportion  to  the  square  root  of  the  Reynolds  number.  This  is  the 
small  scale  that  makes  flows  at  higher  Reynolds  number  more  difficult  to  resolve 
numerically. 

The  Prandtl  number  and  Schmidt  number  for  the  flow  were  chosen  to  be  unity 
(see  section  3.1).  The  ratio  of  specific  heats  was  set  to  7  =  1.4,  and  the  mixture 
fraction,  /,  was  initialized  with  a  hyperbolic  tangent  profile: 

/  =  ^(1  +  tanh(2y)) 

08 


(4.3) 


4.2  Mach  Number  Effects 


In  the  time  developing  mixing  layer  U2  can  be  chosen  as  U2  =  —  l/p2>  so  that 
the  Mach  number  of  each  free-stream  is  the  same  as  the  convective  Mach  number 
(equation  (2.52)).  Thus,  in  the  simulations  we  can  use  free-stream  Mach  number 
and  convective  Mach  number  interchangeably.  Figure  4.5  shows  the  effect  of  Mach 
number  on  the  growth  of  the  most  amplified  disturbance.  The  most  amplified 
wavelength  was  Lx  —  2n/a  —  (7.48,8.06,9.52,13.37)  for  each  of  the  cases  Mi  = 
(0.2, 0.4, 0.6, 0.8).  The  figure  shows  the  growth  in  the  vorticity  thickness  of  the 
layer,  defined  by  equation  (2.21).  For  each  case  a  64  x  81  grid  was  used,  with 
Re  =  400,  Ly  =  10.0  and  disturbance  amplitude  0.05.  Each  wave  grows  until  it 
fills  the  computational  box  and  saturates.  If  the  computation  is  allowed  to  proceed 
further  in  time,  the  behavior  shown  on  figure  4.6,  for  M\  =  0.4,  is  obtained.  The 
layer  thickness  exhibits  a  damped  oscillation  in  time.  Lele  [1989]  showed  that  this 
behavior  is  associated  with  a  ‘nutation’  (shape  change)  of  the  developed  vortex, 
which  also  produces  sound  waves  that  radiate  away  from  the  mixing  layer.  The 
final  structure  is  shown  in  figure  4.7.  The  vorticity  field  is  reminiscent  of  the  Stuart 
vortices  (Stuart  [1967]),  and  the  neutrally  stable  mode  of  linear  analysis  (Michalke 
[1965a]).  The  scalar  is  diffuse,  since  no  new  fluid  is  wrapped  around  the  structure 
in  the  nutation  phase,  only  rotation  of  old  fluid. 

Another  measure  of  thickness  was  tested,  based  on  the  mass-weighted  mean  ve¬ 
locity  profile: 


Figure  4.8  shows  this  measure  compared  with  the  conventional  vorticity  thickness, 
defined  by  equation  4.4.  Two  Mach  numbers  are  considered.  At  M\  =  0.2  there 
is  no  difference  between  the  two  measures,  while  at  M\  =  0.8  the  mass-weighted 
vorticity  thickness  has  a  higher  value  than  the  usual  vorticity  thickness.  It  is  not 
clear  what  the  significance  of  this  difference  is.  In  all  the  following  simulations  the 
usual  vorticity  thickness  is  used. 

The  growth  in  energy  E[kx  =  1)  of  the  disturbances,  computed  by  equation 
(3.53),  is  plotted  on  figure  4.9.  There  is  initially  an  exponential  growth  of  the 


disturbance  energy,  a  straight  line  on  this  semilog  plot.  This  is  followed  by  slower 
growth  with  eventual  saturation  and  decay  in  mode  energy. 

The  effect  of  Mach  number  on  the  developed  structure  is  shown  in  a  series  of 
contour  plots  of  mixture  fraction,  pressure,  vorticity  and  vorticity  divided  by  den¬ 
sity.  These  are  plotted  for  Mach  numbers  M\  =  0.2,  0.4,  0.6  and  0.8  on  figures 
4.10  through  4.13.  All  such  contour  plots  show  equally  spaced  contours  between 
the  maximum  and  minimum  values  shown.  Dashed  contour  lines  are  used  wherever 
a  function  is  negative.  Plots  of  mixture  fraction  are  for  /  —  0.5  and  show  how  fluid 
from  each  of  the  free  streams  is  wrapped  into  the  large-scale  structure.  The  pressure 
fields  show  reduced  pressure  in  the  vortex  cores  (at  z  =  Lz/4,  y  =  0),  and  increased 
pressure  (near  the  isentropic  stagnation  pressure)  at  the  saddle  point  (z  =  3Lz/4, 
y  =  0).  At  M\  —  0.2  the  pressure  is  reduced  by  5.9%  in  the  core,  relative  to  the 
free-stream,  and  raised  by  2.5%  at  the  saddle  point.  At  M\  —  0.8  the  reduction  in 
the  core  is  37.5%,  and  the  rise  at  the  saddle  point  is  48.7%. 

A  comparison  is  made  in  table  4.1  of  the  rise  in  fluid  properties  at  the  saddle 
point,  compared  to  the  rise  assuming  an  isentropic  stagnation  process  (Liepmann 
and  Roshko  [1957])  from  the  free-stream  to  the  stagnation  point: 


T0  ,  ,  7  ~  1  m2 

T,-1+  2  M 1 

(4.5a) 

f  =  (1  + 

Pi  2 

(4.56) 

Po  =  (l  +  7-1M12)1/fr“1> 

Pi  2 

(4.5c) 

where  the  subscript  0  represents  stagnation  conditions.  The  comparison  gets  better 
later  in  the  simulation,  as  the  remaining  vorticity  is  moved  away  from  the  saddle 
point  by  the  local  strain  field.  It  can  be  seen  that,  at  least  for  the  cases  of  equal 
density  considered  here,  the  assumption  of  an  isentropic  process  to  predict  fluid 
conditions  at  the  saddle  point  is  quite  good. 
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Table  4.1  Comparison  of  computed  rise  in  fluid  properties  at  the  saddle 
point  with  rise  assuming  an  isentropic  process. 


Mi 

time 

Ap/Api* 

Ap/Apjg 

Ar/ATjg 

0.2 

15.0 

n 

msm 

1.085 

0.4 

17.6 

Wwm 

1.049 

0.6 

24.0 

.flan 

■am 

1.016 

0.8 

37.2 

0.936 

0.844 

1.170 

The  plots  of  vorticity  and  vorticity  divided  by  density  show  a  clear  change  in 
structure  as  the  Mach  number  is  raised.  The  vortices  become  very  elongated  in  the 
streamwise  direction  at  high  Mach  numbers.  The  effect  is  much  more  noticeable  in 
the  plots  of  vorticity  u>z,  but  is  still  evident  in  the  plots  of  u;2/p.  To  understand  this 
change  in  structure  we  turn  to  the  compressible  vorticity  equations  (2.29)-(2.31). 
Terms  from  these  equations  are  plotted  on  figure  4.14  for  Mi  =  0.6.  The  equation 
for  u)z  has  both  a  dilatational  term  and  a  baroclinic  term  on  the  right  hand  side. 
It  was  found  that  for  mixing  layers  with  equal  free-stream  densities  the  dilatation 
term  was  larger  by  a  factor  of  3  or  4.  This  term  is  plotted  on  figure  4.14c,  and 
shows  a  quadrupole  structure. 

A  physical  explanation  of  the  shape  change  will  now  be  suggested.  A  fluid  element 
approaching  the  structure  from  the  upper  left  hand  side  experiences  an  expanding 
flow  (3uj/3z,-  >  0),  and  a  reduction  in  vorticity,  until  it  is  alongside  the  vortex. 
Then  the  element  is  subject  to  a  compression  (<9u,/dxt  <  0),  with  an  associated 
increase  in  vorticity  as  it  approaches  the  trailing  edge  of  the  vortex,  and  the  stag¬ 
nation  region  behind.  A  similar  process  affects  fluid  elements  approaching  from  the 
lower  right  hand  side.  The  overall  effect  is  to  reduce  vorticity  above  and  below 
the  vortex,  and  to  increase  vorticity  in  front  and  behind  the  vortex,  leading  to  a 
structure  elongated  in  the  streamwise  direction. 

The  same  basic  structure  is  found  for  the  baroclinic  term,  figure  4.l4d,  on  the 
right  hand  side  of  the  equation  for  w2/p.  Hence  this  term  leads  to  elongation  of 
the  w2/p  structure  in  x.  For  this  equation  the  term  on  the  right  hand  side,  though 
large  enough  to  cause  the  elongation  effect,  is  small  enough  to  make  u >2/p  nearly  a 
conserved  quantity.  Fluid  particles  follow  more  closely  the  u>2/p  contours  than  the 
<jjz  contours.  This  is  demonstrated  by  observing  that  the  contours  of  the  mixture 
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fraction,  figure  4.12a,  which  tracks  fluid  elements  during  the  simulation,  show  that 
fluid  has  been  transported  around  paths  resembling  the  u/2/p  contours. 

The  change  in  structure  can  be  related  to  growth  rate  by  considering  the  efficiency 
of  the  various  structures  at  entraining  new  fluid.  Consider  the  limiting  cases  of 
a  circular  vortex,  and  a  fully  elongated  vortex.  The  fully  elongated  vortex  (».e. 
parallel  flow)  does  not  wrap  any  new  fluid  around  the  structure,  and  it  cannot  then 
be  engulfed  and  mixed,  and  there  is  only  growth  by  viscous  diffusion.  The  circular 
vortex  wraps  fluid  from  the  free-streams  around  itself,  and  grows  strongly.  If  we 
assume  a  monotonic  trend  between  these  two  limiting  cases  then  we  see  that  the 
effect  of  an  elongated  vortex  structure  is  to  reduce  the  growth  rate  of  the  mixing 
layer. 


4.3  Effect  of  Mach  Number  on  Pairing 

At  low  Mach  numbers  the  non-linear  mechanism  by  which  the  mixing  layer 
changes  its  lateral  scale  is  observed  in  experiments  to  be  the  vortex  pairing  pro¬ 
cess  (Winant  and  Browand  [1974]).  Two  simulations  were  run,  at  M\  =  0.2  and 
Mi  —  0.6,  with  an  initial  Reynolds  number  of  200,  to  investigate  the  effect  of  Mach 
number  on  pairing.  The  fundamental  ( F )  and  first  subharmonic  (51)  wavelengths 
were  included  in  the  initial  conditions,  with  a  relative  amplitude  of  2:1 

u  =  u  +  A  Real{uf  (y)et(af’I+<^}  +  —  Real{u51(y)e,a5i;r}  (4.6) 

where  <f>  is  the  phase  difference  between  the  two  modes.  Figure  4.15  shows  the 
various  options  for  phase.  The  case  (f>  =  7r/2  was  chosen,  since  this  is  the  optimum  to 
enhance  pairing.  The  opposite  case  is  when  <j>  =  37r/2,  which  corresponds  to  a  slow 
tearing  process,  where  one  vortex  is  trapped  in  the  strain  field  of  the  subharmonic 
and  is  pulled  apart.  All  other  phases  result  in  pairing  or  tearing  in  various  degrees. 

The  time  histories  of  the  vorticity  thickness  are  shown  on  figure  4.16  and  the 
growth  in  energy  on  figure  4.17,  for  Mach  numbers  0.2  and  0.6.  The  energy  plot 
shows  how  the  fundamental  mode  grows  first  and  saturates,  and  then  the  subhar¬ 
monic  wave  takes  over.  The  vorticity  thickness  is  very  sensitive  to  the  structure 
of  the  flow.  Other  measures,  such  as  integrated  momentum  thickness,  or  visual 
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thickness,  do  not  show  all  the  detail  of  the  vorticity  thickness.  Events  in  the  vor- 
ticity  thickness  time  history  can  be  related  to  specific  events  in  the  flow.  This  is 
illustrated  in  the  following  sequence  for  the  M\  =  0.2  simulation. 

(1)  The  fundamental  mode  grows  and  saturates.  Figure  4.18  shows  the  mixture 
fraction,  pressure  and  vorticity  at  time  t  =  9  during  roll-up  of  the  fundamental 
instability. 

(2)  After  saturation  of  the  fundamental  mode,  the  subharmonic  mode  grows  and 
two  of  the  primary  structures  begin  to  rotate  around  each  other.  When  the 
structures  have  rotated  such  that  when  viewed  in  the  x  direction  they  are  two 
separate  vortices,  there  is  an  event  in  the  vorticity  thickness  time  history.  At 
Mi  =  0.2  this  occurs  at  t  =  20,  and  is  illustrated  in  figure  4.19. 

(3)  Rotation  continues  and  the  peak  in  vorticity  thickness  is  reached  when  the 
structures  lie  one  above  the  other.  The  structure  shortly  following  this,  time 
t  =  24  is  shown  on  figure  4.20.  Now  the  structure  on  top  is  beginning  to  move 
downwards,  and  the  two  vortices  are  rotating  around  each  other. 

(4)  The  next  event  in  the  vorticity  thickness  occurs  when  the  vortices  once  again 
lie  above  each  other,  having  rotated  by  180°.  This  occurs  at  time  t  =  27,  and 
is  shown  on  figure  4.21. 

(5)  The  last  point  shown  here  is  at  time  t  =  32,  figure  4.22,  where  the  vortices  have 
rotated  all  the  way  around.  If  the  simulation  is  carried  further,  the  core  of  the 
structure  continues  to  rotate  and  the  vorticity  thickness  continues  a  damped 
oscillatory  behavior,  similar  to  that  found  for  a  single  structure  (figure  4.5).  In 
practice  the  next  subharmonic  would  now  grow  and  the  process  would  repeat, 
beginning  with  step  2  above. 

The  same  sequence  was  observed  at  Mi  =  0.6,  though  the  pairing  process  was 
slowed  both  linearly  and  non-linearly.  Figures  4.23  through  4.26  show  the  structure 
of  the  flow  at  times  16,  29,  37  and  47,  corresponding  approximately  to  items  1  to 
4  in  the  above  pairing  process.  If  we  associate  the  process  between  items  1  and 
2  as  the  linear  growth  of  the  subharmonic,  and  the  remaining  processes  as  non¬ 
linear,  then  it  can  be  seen  from  figure  4.1u  that  both  linear  and  non-linear  aspects 
of  pairing  are  slowed  by  compressibility.  Some  delay  is  expected  due  to  the  increase 
in  wavelength  of  the  most  amplified  disturbances  as  Mach  number  is  increased.  A 
simulation  was  therefore  performed  at  Mi  =  0.2  with  the  same  wavelengths  as  the 
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Ml  =  0.6  case.  Figure  4.27  shows  the  growth  rate  comparison  with  wavelength 
effects  removed,  showing  that  both  linear  and  non-linear  growth  is  indeed  slowed 
by  compressibility. 

The  slowing  of  the  growth  of  the  initial  stage  of  pairing  is  accounted  for  by  the 
effect  of  Mach  number  in  reducing  the  linear  amplification  rate  of  the  subharmonic 
disturbance.  There  are  several  possible  reasons  for  the  reduction  in  growth  rate 
in  the  non-linear  region:  (i)  the  shape  of  the  pairing  vortices  has  changed,  and 
the  more  elongated  structures  may  be  less  inclined  to  pair,  (ii)  the  delay  in  the 
Biot-Savart  type  of  vortex  induction  process,  due  to  the  finite  sound  speed,  may 
make  the  pairing  process  less  efficient,  and  (iii)  continued  compressibility  effects  due 
to  dilatational  and  baroclinic  torques  may  be  slowing  the  rotation  of  the  vortices 
around  each  other.  It  is  not  clear  from  the  current  study  which  of  these  effects  is 
dominant,  or  whether  they  are  all  important. 

The  pairing  involves  a  displacement  of  the  vortices  into  the  free-stream,  which 
distorts  the  streamlines  outside  the  vortical  region,  and  creates  regions  of  locally 
high  velocity.  During  the  pairing  at  =  0.6  it  was  found  that  for  a  short  period  of 
time,  around  the  peak  in  the  vorticity  thickness  curve,  there  was  a  region  of  super¬ 
sonic  flow  above  and  below  the  structures,  although  this  did  not  persist  long  enough 
to  generate  shock  waves.  At  any  higher  Mach  number  shock  waves  would  certainly 
form  during  the  two-dimensional  pairing.  It  is  interesting  to  note  the  evidence  that 
the  flow  would  prefer  not  to  have  shock  waves.  Lele  (private  communication  [1988]) 
has  shown  that  tearing"  may  be  more  common  at  higher  Mach  numbers,  and  Ragab 
and  Wu  [1989]  have  shown  that  a  helical  pairing  may  be  preferred  even  as  low  as 
Mi  =  0.4.  This,  and  the  tendency  toward  three-dimensional  structures  at  high 
Mach  numbers  (see  Chapter  5),  seems  to  indicate  that  the  flow  will  try  to  avoid 
generating  shock  waves  if  at  all  possible.  However,  there  is  no  general  theory  to 
predict  why  this  is  so. 


4.4  Two-Dimensional  Structure 

In  this  section  the  two-dimensional  structure,  under  the  influence  of  compress¬ 
ibility,  is  explored.  The  evolution  of  the  vorticity  and  scalar  fields  at  this  Mach 
number  have  already  been  described  in  section  4.2.  Here  we  consider  the  behavior 
of  other  fluid  properties,  such  as  stagnation  enthalpy,  entropy,  and  strain  rate. 
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Contours  of  temperature  at  M\  =  0.6  are  shown  on  figure  4.28a.  They  show 
high  temperature  at  the  stagnation  point  in  the  braid,  and  lower  temperature  in 
the  vortex  core.  The  stagnation  enthalpy  for  a  perfect  gas  with  constant  specific 
heats  is  given  by: 


and  non-dimensionally: 


*2  *2 

*  *  *  U  4-  V 

h'  =  cjr*  +  — J— 


H  = 


C2  U2  +  V2 


(7-1) 


+ 


(4.7a) 


(4.76) 


A  contour  plot  of  H  is  shown  on  figure  4.286.  The  stagnation  enthalpy  is  lowest  in 
the  vortex  core,  and  remains  approximately  constant  in  the  rest  of  the  flow.  This 
effect  can  be  understood  by  considering  an  inviscid,  non-conducting  form  of  the 
energy  equation  (Liepmann  and  Roshko  [1957]). 


DH  1  dp 
Dt  p  dt 

In  the  cores  of  the  structures  during  roll-up  the  pressure  is  reducing  with  time,  so  a 
fluid  particle  being  wrapped  around  the  developing  structure  experiences  decreasing 
pressure,  and  by  equation  (4.8)  the  stagnation  enthalpy  drops.  In  situations  where 
the  two  free  streams  have  different  stagnation  enthalpies  the  reduced  H  in  the  core 
is  overshadowed  by  the  mean  H  profile. 

The  entropy  of  the  flow  (non-dimensionalized  as  s  =  s*T?  /Uf  since  e  is  normal- 
2 

ized  by  ),  relative  to  the  free-stream  1,  is  defined  by 


_  In  T  In  p 

7(7  -  1  )Af2  7 M2 

A  contour  plot  is  shown  on  figure  4.28c.  The  entropy  structure  is  similar  to  the 
vorticity,  with  high  entropy  found  in  the  center  of  the  vortex.  The  assumption  of 
an  isentropic  process  from  the  free  streams  to  the  saddle  point  can  be  seen  to  be 
valid  for  free-streams  with  equal  entropy. 
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Strain  rate  in  the  flow  is  of  interest  from  two  perspectives.  Firstly,  because  it 
strongly  affects  the  development  of  streamwise  vorticity  in  the  flow,  and  strain  rate 
near  the  saddle  point  can  lead  to  the  ‘collapse’  of  streamwise  vorticity  into  stream- 
wise  vortices  (Lin  and  Corcos  [1984]).  This  produces  mushroom  structures  in  the 
scalar  field,  observed  experimentally  by  Bernal  and  Roshko  [1986].  Secondly,  strain 
rate  is  of  interest  when  applications  to  combustion  are  considered.  If  the  strain  rate 
is  too  high  in  any  region  of  the  flow  then  any  flame  that  forms  may  be  quenched, 
and  reactions  generating  heat  release  may  be  prevented.  Unlike  incompressible  flow, 
the  strain  rate  tensor  is  not  trace-free  and  there  are  alternative  ways  to  document 
the  strain  rate  magnitude.  The  usual  definition  of  strain  rate  tensor  is 

=  +  (4-10) 

An  anisotropic  strain  rate  tensor  can  also  be  defined  (Reynolds  [1988]): 

c-a  _  c  Skk&ij 

-  bi]  3 

The  magnitudes  of  the  strain  rates  for  each  case  are  defined  by: 

|S‘I  =  v/spS  (4-13) 

For  the  simulations  presented  here  it  was  found  that  the  high  strain  rate  regions 
in  the  flow  were  not  associated  with  regions  of  high  dilatation,  so  the  values  for 
peak  strain  rate  magnitude  were  not  affected  by  choice  of  equation  (4.12)  or  (4.13). 
At  different  times  in  the  simulation  for  Mi  =  0.6  contour  plots  of  the  strain  rate 
field  are  shown  on  figures  4.29  and  4.30,  and  slices  through  the  flow  at  various  x 
locations  are  shown  on  figure  4.31.  It  can  be  seen  that  the  usual  picture  of  strain 
rate,  high  in  the  braid  regions  and  low  in  the  cores,  is  observed  at  time  t  =  18.2 
during  the  early  stage  of  roll-up.  However  at  later  times,  e.g.  at  t  =  24.0,  figure 
4.316,  the  strain  rate  peak  in  the  braid  is  not  dominant,  and  regions  of  high  strain 
rate  are  also  found  within  the  vortex  core. 

For  combustion  application,  it  would  be  useful  to  be  able  to  model  the  magnitude 
of  the  peak  strain  rate  in  the  flow.  The  time  history  of  the  peak  strain  rate  is  shown 


(4.11) 

(4.12) 
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on  figure  4.32.  When  normalized  by  the  initial  vorticity  thickness  and  free-stream 
velocity,  the  strain  rate  is  seen  to  decrease  during  the  simulation.  When  the  strain 
rate  is  renormalized  by  the  local  vorticity  thickness  and  velocity  difference,  shown 
on  curve  2  of  figure  4.32,  we  find  that  the  strain  rate  increases  slightly  during  roll¬ 
up,  and  then  declines.  In  the  latter  case  the  non-dimensional  peak  strain  rate  is 
always  of  order  unity,  so  that  peak  strain  rate  in  the  mixing  layer  can  be  modeled 
by  the  ratio  of  the  free-stream  velocity  difference  to  the  local  vorticity  thickness. 

4.5  Effect  of  Density  Ratio 

Most  compressible  mixing  layer  experiments  are  performed  for  free-streams  with 
unequal  densities.  Additional  baroclinic  torques  can  appear  in  the  flow,  especially 
in  the  braid  region  where  the  two  free-streams  are  brought  close  together  and  steep 
density  gradients  are  formed.  To  investigate  these  effects,  two  simulations  are  com¬ 
pared;  one  has  free-streams  of  equal  density,  and  the  other  has  a  density  ratio 
P2  =  0.5  (T2  =  2). 

The  simulations  are  run  at  Mc  =  0.6,  with  a  Reynolds  number  of  400.  To  check 
the  convective  velocity  formula,  we  arrange  U2  so  that  the  predicted  convective 
velocity  is  zero.  For  the  density  ratio  of  0.5,  this  means  setting  U2  =  —  y/2.  The 
velocity  profile  is  then  obtained  from  the  error  function  as: 

u  =  ^  i  +  U2  +  (1  -  ^2)erf(yV7r)]  (4.14) 

The  temperature  profile  is  obtained  from  the  velocity  profile  in  the  usual  way  (equa¬ 
tion  (2.13)). 

The  growth  in  vorticity  thickness  for  each  simulation  is  shown  on  figure  4.33, 
and  the  growth  in  integrated  energy  £(1)  (the  fundamental  mode),  on  figure  4.34. 
The  plots  of  E(l)  show  that  the  simulations  start  from  a  different  origin.  This 
is  simply  a  reflection  of  differences  in  the  shape  of  the  eigenfunctions,  since  the 
normalization  is  done  by  amplitude  of  the  u  eigenfunction,  and  not  by  integrated 
energy.  The  developed  structure  is  shown  on  figure  4.35,  and  can  be  compared  with 
the  simulation  for  uniform  density  shown  on  figure  4.12.  It  can  be  seen  that  there  is 
generation  of  vorticity  with  the  opposite  sign  to  that  of  the  dominant  roll-up,  which 
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can  be  explained  by  considerations  of  the  baroclinic  torque  acting  in  the  saddle- 
point  region  of  the  mixing  layer.  A  sketch  of  the  saddle  region  is  shown  on  figure 
4.36.  In  region  1  the  gradients  of  pressure  and  density  are  nearly  perpendicular 
to  each  other  and  application  of  a  right  hand  rule  shows  that  clockwise  vorticity 
(positive  in  the  current  definition)  is  produced.  This  is  then  convected  away  from 
the  saddle  point  by  the  local  strain  rate  field,  and  becomes  wrapped  around  the 
structure  to  the  left.  In  region  2  the  pressure  and  density  gradients  act  to  produce 
vorticity  with  the  opposite  sign  to  the  main  roll-up  (anti-clockwise),  which  becomes 
wrapped  around  the  structure  to  the  right.  The  counter-clockwise  vorticity  does  not 
produce  vortices  of  the  opposite  sign,  it  merely  modifies  the  path  of  fluid  particles 
rotating  around  the  main  structure,  and  results  in  a  more  asymmetric  scalar  field. 

A  simulation  was  run  with  a  density  ratio  of  0.2  to  investigate  the  Uc  formula 
(1.2),  compared  with  the  predictions  of  the  linear  instability  theory.  For  this  case 
(T2  =  5,  U2  =  —  \/5)  the  convective  velocity  formula  predicts  zero  convective  veloc¬ 
ity,  whilst  from  the  linear  theory  the  phase  speed  of  the  growing  wave  is  0.05  and 
the  phase  speed  of  the  neutral  mode  is  0.16,  which  are  both  very  definitely  non-zero. 
The  developed  structure  is  shown  in  figures  4.37.  The  scalar  field  is  highly  asym¬ 
metric  due  to  the  baroclinic  terms.  Figure  4.38  shows  the  initial  and  final  pressure 
contours.  Initially  the  pressure  minima  is  at  x  =  Lx/A  and  the  pressure  maxima  is 
at  x  =  3Lz/4.  At  the  later  time  the  pressure  minima  has  moved  to  x  =  Lx/ 2  and 
the  pressure  maxima  is  at  x  =  Lx/§.  The  structure  has  moved  to  the  right.  The 
approximate  convective  velocity,  based  on  the  movement  of  the  stagnation  point 
(a  peak  in  presage)  is  0.15.  It  is  0.12  based  on  the  movement  of  the  minima  of 
pressure  in  the  vortex  core.  The  reason  for  the  difference  is  that  the  structure  is 
evolving  and  the  convective  velocity  changes  depending  on  which  feature  of  the  flow 
we  look  at.  These  rough  estimates  agree  better  with  the  phase  speed  of  the  neutral 
mode  than  with  the  LJC  formula  (1.2),  which  would  predict  no  convection. 


4.6  Embedded  Shock  Waves 

In  two-dimensional  simulations  above  a  convective  Mach  number  of  0.7  it  was 
found  that  shock  waves  developed  in  the  flow.  Pressure  contours  for  a  simulation 
at  Mj  =  0.8  at  a  Reynolds  number  of  400  are  shown  on  figure  4.39a,  and  show 
the  location  of  the  shock  waves.  Local  Mach  contours  are  shown  on  figure  4.396. 
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Flow  around  the  top  of  the  vortex  is  accelerated  to  supersonic  speeds,  and  then  has 
to  stagnate  at  the  saddle  point  in  the  braid.  It  does  this  by  compression  through 
a  weak  (Mmax  «  1.2)  near-normal  shock  wave.  The  process  is  similar  to  the  flow 
around  a  transonic  airfoil.  The  same  process  occurs  for  the  lower  stream  moving 
from  right  to  left  below  the  vortex.  Profiles  of  pressure,  density  and  temperature 
through  the  shock  wave  are  shown  on  figure  4.40.  They  show  that  for  this  flow 
condition  the  shock  wave  is  adequately  captured  by  the  numerical  method,  using  a 
grid  with  192  points  in  the  streamwise  direction. 

The  computational  problem  is  that  at  lower  Reynolds  numbers  ( t.g .  200)  the 
roll-up  of  the  vortex  is  not  strong  enough  to  generate  a  shock  wave  (i.e.  we  are  in  a 
viscous  regime),  while  at  higher  Reynolds  numbers  the  shock  is  stronger,  and  thin¬ 
ner,  so  that  a  prohibitively  large  number  of  modes  are  required  to  resolve  it  properly. 
During  an  earlier  phase  of  this  project,  TVD  (Total  Variation  Diminishing)  numeri¬ 
cal  methods  were  tested  for  mixing  layer  flows.  These  numerical  methods  are  of  the 
shock-capturing  genre  (Yee  [1989]),  and  were  found  to  do  a  good  job  of  capturing 
the  shock  waves,  no  matter  how  strong,  without  encountering  numerical  difficulties. 
However,  this  shock  capturing  capability  is  obtained  by  introducing  extra  dissipa¬ 
tion  in  regions  of  the  flow  where  oscillations  develop.  This  extra  dissipation  also 
tends  to  damp  out  the  growth  of  the  large-scale  structure  in  the  flow.  It  was  found 
that  with  the  TVD  schemes  a  large  number  of  grid  points  were  required  in  order  to 
capture  the  instability  correctly  (Sandham  and  Yee  [1989]).  Neither  spectral  meth¬ 
ods  (relying  on  shock  resolution) ,  nor  TVD  methods  (shock  capturing)  were  found 
to  be  efficient  for  flows  involving  both  shock  waves  and  growth  of  instabilities. 

A  model  two-dimensional  large-scale  structure  with  shock  waves  has  been  pro¬ 
posed  by  Papamoschou  [1989],  and  separately  by  Dimotakis  [1989],  to  explain  some 
of  Papamoschou’s  experimental  results  concerning  convective  velocities.  (An  al¬ 
ternative  explanation  of  the  results  using  linear  stability  theory  was  presented  in 
section  2.4.)  The  shock  waves  proposed  by  Papamoschou  are  strong  and  are  in¬ 
clined  to  the  streamwise  direction,  so  that  there  is  a  large  entropy  rise  through 
the  shock  wave.  This  breaks  the  symmetry  in  the  convective  velocity  formula,  and 
was  offered  as  an  explanation  of  the  experimental  finding  that  convective  speeds 
are  biassed  towards  one  free-stream  speed  or  the  othrr  (Papamoschou  [1989]).  The 
model  structure  made  no  prediction  of  the  angle  or  strength  of  the  shock  waves. 
However,  the  shocks  observed  in  the  current  work  were  always  weak,  nearly  normal 
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to  the  streamwise  direction,  and  were  found  on  both  sides  of  the  vortices.  Thus, 
they  are  not  reminiscent  of  the  model  shock  structure  proposed  by  Papamoschou 
and  by  Dimotakis. 


4.7  Simulation  of  Supersonic  Mode  Instability 

The  ‘radiating  vorticity’  modes  of  instability,  supersonic  with  respect  to  one  of 
the  free-streams,  are  the  most  amplified  modes  in  two  dimensions  above  a  convective 
Mach  number  of  approximately  1  (see  section  2.2.4).  Although  we  expect  three- 
dimensional  modes  to  be  dominant  at  these  Mach  numbers,  it  is  interesting  to 
simulate  the  non-linear  development  of  these  modes,  which  might  be  excited  in 
experiments  by  strong  two-dimensional  forcing  of  the  mixing  layer. 

The  time-developing  mixing  layer  with  equal  densities  at  a  Mach  number  of  1.05 
was  chosen  for  simulation.  The  Reynolds  number  was  fixed  at  800,  which  is  just 
high  enough  so  that  the  inviscid  eigenfunction  is  amplified.  The  high  Reynolds 
number  necessary  for  simulation  of  this  weak  instability  means  that  shock  waves 
cannot  resolved,  since  these  will  be  very  thin.  The  computations  have  to  be  stopped 
as  soon  as  shock  waves  form. 

The  mode  that  is  supersonic  with  respect  to  the  lower  stream  was  simulated.  The 
other  mode,  supersonic  relative  to  the  upper  stream,  will  have  the  same  behavior. 
From  the  linear  stability  theory  the  amplification  rate  of  this  mode  is  29%  of  the 
amplification  rate  of  the  most  amplified  three-dimensional  wave  (angle  56°). 

The  growth  in  vorticity  thickness  and  energy  E(  1)  for  this  supersonic  mode  are 
shown  in  figures  4.41  and  4.42  respectively.  Contours  of  relevant  features  are  shown 
on  figure  4.43  at  time  t  =  82.  We  see  that  roll-up  does  occur,  so  that  this  is  indeed 
a  ‘vortical’  instability  mode,  which  can  act  to  mix  fluid  from  the  two  streams. 
However,  the  roll-up  is  very  weak  and  takes  a  long  time  to  develop.  The  structure, 
shown  for  example  in  the  contours  of  density  on  figure  4.436,  seems  to  reside  on  the 
upper  side  of  the  mixing  layer,  and  the  contours  of  mixture  fraction  show  that  the 
cores  of  the  structures  are  mostly  fluid  from  the  upper  stream,  i.e.  the  stream  that  is 
subsonic  relative  to  the  structure.  The  pressure  contours,  figure  4.43c,  show  clearly 
the  radiating  nature  of  the  instability.  On  the  upper  side  we  have  the  usual  picture 
of  pressure:  high  at  the  stagnation  point  and  low  in  the  vortical  structure.  However 
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on  the  lower  side,  where  the  free-stream  is  supersonic  relative  to  the  structure, 
we  find  an  arrangement  of  expansion  and  compression  waves.  The  compression 
waves  develop  into  shock  waves  during  the  simulation.  At  a  later  time  t  =  92.5 
the  structure,  shown  on  figure  4.44,  has  developed  shock  waves  which  can  no  longer 
be  resolved  on  the  computational  grid;  oscillations  can  be  observed  in  the  pressure 
contours  in  the  lower  left  corner  of  the  plot. 

The  presence  of  walls  in  experiments  would  cause  the  waves  to  reflect  back  into 
the  shear  layer.  Such  a  mechanism  lies  at  the  heart  of  the  wall-modes  and  Mack 
modes  proposed  by  other  authors  (Greenough  et  a/.[l989],  Mack  [1989]),  though  is 
is  not  clear  that  those  modes  are  vortical  modes  that  could  act  to  mix  fluid.  In  any 
case,  as  long  as  the  walls  aren’t  too  close,  the  three-dimensional  instability  is  more 
amplified,  and  this  is  the  subject  of  the  next  chapter. 


4.8  Chapter  Summary 

The  main  effect  of  compressibility  in  two-dimensional  simulations  is  to  damp 
the  growth  of  the  instability,  both  linearly  and  non-linearly.  As  Mach  number 
is  increased  the  developed  vortical  structure  becomes  elongated  in  the  streamwise 
direction  and  less  efficient  at  wrapping  new  fluid  into  the  mixing  layer.  The  pairing 
process  is  also  slowed  by  compressibility.  For  convective  Mach  numbers  above  0.7 
we  find  weak  shock  waves  developing  in  the  flow. 

Above  a  convective  Mach  number  of  1  the  only  unstable  two-dimensional  modes 
in  the  mixing  layer  are  the  weakly  amplified  supersonic  modes.  Simulation  of  these 
modes  showed  that  the  structures  developing  from  these  instabilities  are  weakly 
vortical,  and  are  able  to  mix  fluid.  The  structures  are  supersonic  relative  the  one 
of  the  free-streams,  and  on  that  side  a  pattern  of  shocks  and  expansion  fans  forms. 
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CHAPTER  5 

Three-Dimensional  Simulations 

In  this  chapter  simulations  are  presented  of  the  three-dimensional  instabilities, 
which  are  expected  to  dominate  the  mixing  layer  at  high  Mach  number.  The  large- 
scale  structures  which  develop  from  the  non-linear  growth  of  these  instabilities  are 
analyzed. 

5.1  Initial  Conditions  and  Parameters 


The  time-developing  mixing  layer  is  again  chosen  for  simulation.  The  initial 
mean  flow  is  the  same  as  for  the  two-dimensional  simulations  of  Chapter  4,  with 
the  mean  velocity  given  by  equation  (4.1)  and  the  mean  temperature  by  equation 
(2.13). 

Three  Mach  numbers  were  selected  for  detailed  study.  A  free-stream  Mach 
number  M\  =  0.4  was  chosen  as  a  low  Mach  number  where  we  expect  to  find 
near-incompressible  behavior.  An  intermediate  Mach  number  was  chosen  to  be 
M\  =  0.8.  At  this  Mach  number  an  oblique  wave  is  slightly  more  amplified  than 
the  two-dimensional  wave,  and  a  broad  range  of  waves  in  between  are  about  equally 
amplified.  Mach  number  M\  —  1.05  was  selected  as  the  high  Mach  number  case  for 
study. 

The  instability  characteristics  of  the  flow  change  around  a  convective  Mach  num¬ 
ber  of  1.  Just  below  Mc  =  1  the  two-dimensional  wave,  though  no  longer  the  most 
amplified  wave,  is  still  only  a  factor  of  about  2  less  amplified  than  the  most  ampli¬ 
fied  wave,  and  may  still  be  expected  to  ,  \ ,  an  effect  on  the  development  of  the 
flow.  Above  Mc  =  1  the  two-dimensional  w  ,  is  much  more  weakly  amplified  than 
the  oblique  waves,  and  the  flow  is  expected  to  be  dominated  by  a  narrow  band  of 
oblique  waves.  This  change  is  illustrated  on  figure  5.1,  where  the  linear  amplifica¬ 
tion  rate  is  plotted  as  a  function  of  angle  at  Mach  numbers  Mj  =  0.95  and  1.05  for 
the  time-developing  mixing  layer  with  equal  free-stream  densities.  The  simulation 
at  M\  =  1.05  is  expected  to  be  typical  of  higher  Mach  number  mixing  layers,  since 


the  only  further  change  in  the  stability  characteristics  of  the  flow  is  the  increasing 
obliquity  of  the  most  amplified  disturbances. 

The  Reynolds  numbers  for  the  simulations  were  chosen  by  considering  the  effect 
of  viscosity  on  the  growth  rate  of  an  eigenfunction  from  inviscid  linear  stability 
analysis.  Figure  5.2  shows  the  effect  of  Reynolds  number  on  the  amplification  rate 
for  three  different  oblique  waves:  a  45°  wave  at  M\  =  0.4,  a  45°  wave  at  M\  =  0.8, 
and  a  60°  wave  at  M\  =  1.05.  The  latter  two  are  expected  to  be  the  most  amplified 
waves  at  their  respective  Mach  numbers.  Ideally,  we  would  like  to  be  at  the  high 
Reynolds  number  end  of  these  curves  in  order  to  capture  the  inviscid  nature  of 
the  instabilities.  The  Reynolds  number  was  selected  to  be  400,  600  or  800  for  the 
Mach  numbers  0.4,  0.8  and  1.05  respectively.  The  problem  with  simulating  high 
Mach  number  flows  is  the  higher  Reynolds  numbers  that  are  required  to  capture 
the  instability.  The  higher  the  Reynolds  number  the  more  modes  are  required  later 
in  the  simulation  in  order  to  resolve  the  flow.  For  this  reason  the  highest  Mach 
number  simulated  was  1.05. 

The  three-dimensional  simulations  presented  in  this  chapter  were  all  run  for 
a  Prandtl  number  of  1  and  a  Schmidt  number  of  1.  Initial  perturbations  were 
usually  specified  as  eigenfunctions  from  the  linear  stability  analysis.  However,  some 
simulations  (sections  5.4  and  5.6)  were  performed  with  random  initial  conditions,  to 
check  that  the  most  amplified  waves  were  indeed  those  predicted  by  linear  stability 
analysis.  The  simulations  were  usually  started  on  a  16  x  99  x  16  grid,  with  boundary 
conditions  applied  at  y  =  ±5  initial  vorticity  thicknesses  away  from  the  mixing  layer 
centerline.  As  the  high  wavenumbers  gained  energy  the  number  of  points  used  in  x 
and  z  was  extended,  ending  at  typically  96  x  99  x  96. 


5.2  Simulations  at  Low  Mach  Number 

At  Mi  =  0.4  we  expect  to  get  nearly  incompressible  behavior.  This  case  there¬ 
fore  partially  serves  as  a  check  on  the  code,  comparing  against  previous  incompress¬ 
ible  mixing  layer  simulations  and  against  incompressible  secondary  stability  theory. 
However,  some  compressibility  effects  may  be  evident  even  at  this  low  Mach  num¬ 
ber  since  Ragab  and  Wu  [1989]  found  some  changes  from  incompressible  secondary 
stability  behavior  at  M\  =  0.4. 
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The  major  work  in  secondary  stability  analysis  for  the  incompressible  mixing 
layer  is  that  of  Pierrehumbert  and  Widnall  [1982].  They  used  Stuart  [1967]  vortices 
as  the  base  flow  and  solved  the  resulting  Floquet  problem  for  both  fundamental  and 
subharmonic  instabilities.  The  fundamental  instabilities  were  of  the  core  ‘bulging’ 
or  core  ‘translative’  kinds,  while  the  subharmonic  instabilities  were  either  regular 
pairing  or  a  ‘helical  pairing’.  Simulations  of  all  these  instabilities  can  be  made  by 
choosing  an  initial  disturbance  field  of  a  two-dimensional  wave,  and  two  equal  and 
opposite  oblique  waves. 

The  fundamental  mode  instabilities  have  the  same  streamwise  wavelength  as  the 
initial  roll-up.  The  initial  disturbance  field  is  specified  by: 

u'  =  ^iReal{fi(a,0)e*l“‘+*>}  +  -42Rea]{S(a,£)e*(“I+W  +  “(«. 

(5.1) 

where  u(a,/?)  is  an  eigenfunction  of  the  linear  instability  wave  with  streamwise 
wavelength  Lx  =  2 n/a  and  spanwise  wavelength  Lz  =  27r//?.  Similar  disturbances 
are  added  for  p',  v',  w'  and  T' .  The  phase  of  the  oblique  waves  relative  to  each 
other  is  not  important,  since  this  does  not  change  the  basic  pattern  of  the  addition 
of  two  oblique  waves,  only  translating  it  in  space.  However  the  phase,  <£,  of  the 
two-dimensional  wave  relative  to  the  pair  of  oblique  waves  is  important. 

With  the  phase  <f>  —  0  we  get  the  ‘bulging’  mode  of  Pierrehumbert  and  Widnall 
[1982].  The  structure  of  the  initial  field  for  this  case  can  be  shown  by  a  cut  through 
the  x-z  plane  at  y  =  0,  figure  5.3.  Pressure  contours  are  shown  on  figure  5.3 a 
and  contours  of  u>x,  u>y  and  u>z  on  figures  5.3 b-d  respectively.  The  pressure  minima, 
which  shows  where  the  core  of  the  vortical  structure  will  form,  is  at  x  =  Lx/ 4,  and  is 
fairly  uniform  across  the  span.  The  contours  of  u>x  and  u>y  show  a  pattern  of  counter¬ 
rotating  vortices,  resulting  from  the  addition  of  the  two  equal  and  opposite  oblique 
waves.  For  these  plots  the  sign  of  vorticity  has  been  chosen  so  that  clockwise  motion 
is  shown  with  solid  contours  and  counter  clockwise  motion  with  dashed  contours. 
It  can  be  seen  that  the  effect  of  these  vortices  is  to  distort  the  core  of  the  two- 
dimensional  instability  and  give  it  a  varying  cross-sectional  area  in  the  spanwise 
direction. 

With  the  phase  4>  changed  to  x/2  we  get  the  ‘translative’  mode  of  Pierrehumbert 
and  Widnall  [1982].  The  initial  condition  is  shown  on  figure  5.4.  Now  the  primary 
roller  is  at  x  =  0,  which  coincides  with  a  peak  in  both  streamwise  and  vertical 
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components  of  vorticity.  The  effect  is  to  make  the  pressure  minima  oscillate  in  the 
spanwise  direction,  alternately  forcing  it  forwards  or  backwards. 

Direct  simulations  were  made  for  both  the  bulging  and  translative  modes  at 
Mi  =  0.4.  The  Reynolds  number  was  set  to  400  and  the  disturbance  amplitudes 
were  A\  =  0.05  and  A<i  =  0.025  (equation  (5.1)).  The  angle  for  the  most  amplified 
disturbance  from  Pierrehumbert  and  Widnall  was  56.3°,  but  the  variation  with 
angle  was  weak.  It  was  found  that  at  this  angle  and  Reynolds  number  the  initial 
instability  wave  was  damped.  A  smaller  angle  is  needed  in  order  to  get  a  disturbance 
that  grows  on  the  mean  flow  and  45°  was  selected.  Plots  of  the  vorticity  thickness 
growth  for  the  two  simulations  are  shown  on  figure  5.5,  and  the  time  history  of  the 
energy  E  (equation  (3.53))  is  shown  on  figure  5.6  for  the  bulging  mode  and  on  figure 
5.7  for  the  translative  mode.  Modes  are  defined  by  (kx,kz)  where  kx  and  kz  are 
integer  wavenumbers  in  the  x  and  z  directions.  Mode  (1,0)  is  the  two-dimensional 
wave,  and  (1,1)  and  (1,-1)  are  the  two  oblique  waves.  For  these  simulations  the 
oblique  waves  are  equal  and  opposite,  and  grow  in  exactly  the  same  way.  There 
is  a  line  of  symmetry  in  the  simulations  at  z  =  Lz/ 2,  the  preservation  of  which 
serves  as  another  check  on  the  code.  Comparison  of  figures  5.6  and  5.7  shows  that 
the  (1,0)  mode  is  unaffected  by  its  phasing  relative  to  the  oblique  waves.  However 
the  growth  of  the  three-dimensional  waves  is  much  stronger  for  the  phase  <f>  —  n/2 
than  for  <j>  =  0,  which  agrees  with  Pierrehumbert  and  Widnall’s  finding  that  the 
translative  mode  of  instability  is  more  amplified  than  the  bulging  mode. 

The  best  method  that  was  found  for  identification  of  rotational  regions  in  the 
flow  was  to  search  for  a  local  minima  in  pressure.  Perspective  views  of  a  surface  of 
constant  pressure,  enclosing  a  region  of  low  pressure,  are  shown  on  figures  5.8  and 
5.9  for  the  bulging  and  translative  modes  respectively.  The  bulging  mode  appears  to 
be  only  weakly  unstable,  and  the  developed  structure  is  very  two-dimensional.  The 
translative  mode  results  in  a  vortex  core  that  oscillates  in  the  spanwise  direction. 

The  fate  of  the  streamwise  and  vertical  vorticity  initially  in  the  saddle  point 
region  of  high  strain  between  two  large  spanwise  rollers  is  of  much  interest,  since 
Lin  and  Corcos  [1984]  showed  that  straining  of  vorticity  in  this  region  can  lead  to  the 
formation  of  streamwise  vortices.  These  vortices  produce  mushroom-like  structures 
in  a  scalar  field,  as  observed  in  experiments  by  Bernal  and  Roshko  [1986].  For  the 
initial  condition  of  the  bulging  mode  (figure  5.3)  the  saddle  point  lies  between  two 
of  the  regions  of  streamwise  vorticity.  Development  of  the  primary  roller  creates 
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a  straining  in  this  region,  which  pulls  the  streamwise  and  vertical  vorticity  away 
from  the  saddle  point.  A  perspective  view  of  the  developed  streamwise  vorticity 
is  shown  on  figure  5.10,  corresponding  to  the  pressure  field  of  figure  5.8.  There  is 
no  streamwise  vorticity  in  the  saddle  region,  and  the  streamwise  vorticity  from  the 
initial  condition  has  been  pulled  away  from  the  saddle  point  towards  the  vortex 
core. 

In  the  initial  condition  for  the  translative  mode  (figure  5.4)  there  is  streamwise 
and  vertical  vorticity  centered  on  the  saddle  point  (located  by  the  pressure  peak  at 
z  =  Lx/2).  Now  the  straining  field  acts  to  pull  this  vorticity  into  long  thin  regions  of 
vorticity  (the  braids).  If  the  initial  vorticity  is  strong  enough  Lin  and  Corcos  [1984] 
predict  a  ‘collapse’  of  the  streamwise  and  vertical  vorticity  into  vortices  aligned  with 
the  principal  axis  of  strain.  This  effect  has  been  confirmed  in  the  incompressible 
numerical  simulations  of  Rogers  and  Moser  [1989].  Figure  5.11  shows  the  streamwise 
vorticity  for  the  structure  that  develops  from  the  translative  instability  at  M\  = 
0.4.  The  vorticity  has  indeed  become  elongated  in  the  saddle  region,  but  for  this 
initial  condition  there  has  been  no  collapse  into  the  near-circular  vortices  of  Lin  and 
Corcos.  A  higher  amplitude  of  initial  disturbances  would  be  required  in  order  to 
get  these  structures. 

The  subharmonic  instabilities  of  Pierrehumbert  and  Widnall  [1982]  are  those 
with  wavelength  twice  that  of  the  primary  roll-up.  The  two-dimensional  subhar¬ 
monic  instability  is  the  two-dimensional  pairing  process,  described  in  section  4.3 
of  this  work.  The  other  subharmonic  instability  was  labeled  a  ‘helical  pairing’  by 
Pierrehumbert  and  Widnall,  and  has  not  been  simulated  numerically  before.  This 
mode  can  be  obtained  by  the  following  combination  of  instability  waves: 

v!  =  AiReal{u(2o,0)e*<2aI+*>}  +  y42Real{u(a,/J)ti(“+W  + 

(5.2) 

The  computational  box  length  is  Lx  =  27r/a  and  the  box  width  is  Lz  =  2? r//?.  The 
initial  field  with  phase  <f>  =  x/2  is  shown  on  figure  5.12.  The  spanwise  vortices  are 
located  by  the  pressure  minima  at  x  =  0,  Lx/2.  The  ‘helical’  instability  occurs  when 
the  streamwise  and  vertical  vorticity  from  the  oblique  waves  is  superimposed  on  the 
centers  of  these  structures.  At  z  =  Lz/2  the  left  structure  is  lifted  and  moved  to 
the  right,  while  the  right  structure  is  pushed  down  and  moved  to  the  left.  At  z  =  0 
the  situation  is  reversed. 
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A  ‘helical’  pairing  simulation  was  run  with  an  initial  Reynolds  number  of  200,  a 
box  size  of  Lx  =  Lz  =  15.7  (».e.  the  45°  oblique  subharmonic  wave  was  excited),  and 
boundaries  at  y  =  ±10.  The  growth  in  vorticity  thickness  is  shown  on  figure  5.13 
and  the  growth  in  mode  energies  on  figure  5.14.  Mode  (2,0)  is  the  two-dimensional 
fundamental  wave,  and  modes  (1,1)  and  (1,-1)  are  the  oblique  subharmonic  waves. 
The  structure  after  roll-up  of  the  primary  instability  is  shown  as  a  perspective  view 
of  pressure  on  figure  5.15.  The  two  vortices  have  been  perturbed  by  the  action  of 
the  oblique  waves.  At  a  later  time  the  structure  is  shown  on  figure  5.16.  The  per¬ 
turbation  shape  has  developed  into  strong  oscillation  of  the  vortices  in  the  spanwise 
direction.  However,  contrary  to  Pierrehumbert  and  Widnall’s  interpretation,  we 
find  no  evidence  for  rotation  of  the  two  vortices  around  each  other,  in  other  words 
no  merging  by  pairing.  Instead  the  flow  has  evolved  into  a  subharmonic  structure, 
with  wavelength  twice  that  of  the  original  roll-up.  The  final  structure  resembles  the 
hairpin  structures  found  in  wall  boundary  layer  transition,  though  here  there  is  an 
antisymmetry  between  the  upper  and  lower  parts  of  the  vortex  tubes,  and  there  are 
4  heads  and  4  legs  per  periodic  structure.  Near  the  head  of  the  hairpins  there  is  an 
induced  motion  due  to  the  legs,  which  will  tend  to  pull  the  structure  into  a  more 
upright  orientation,  and  will  oppose  any  tendency  to  rotate  around  neighboring 
structures.  This  prevents  the  ‘helical’  instability  becoming  a  ‘helical  pairing’.  The 
presence  of  any  two-dimensional  wave  would  eventually  lead  to  pairing,  but  if  the 
time  to  reach  this  pairing  is  long  compared  to  the  time  taken  for  the  next  stage  of 
instability  to  develop  we  may  see  structures  similar  to  those  on  figure  5.16. 

Cuts  through  the  mixture  fraction  and  pressure  fields  at  the  plane  y  =  0  are 
shown  on  figure  5.17.  The  final  structure  found  here  is  very  similar  to  that  which 
will  be  presented  for  high  Mach  number  flows  in  section  5.5.  The  difference  is  that 
here  the  oblique  waves  grow  on  a  base  flow  of  developed  spanwise  rollers,  whereas 
in  section  5.5  the  oblique  waves  grow  on  the  unperturbed  mean  flow. 

Evidence  for  the  existence  of  these  structures  in  experiments  is  slim,  but  a  ‘he¬ 
lical’  structure  was  claimed  to  have  been  observed  by  Chandrsuda  et  al.  [1978], 
which  may  have  been  due  to  an  underlying  structure  like  the  one  in  figure  5.16. 
For  incompressible  flows  this  instability  is  less  amplified  than  the  classical  two- 
dimensional  pairing,  so  it  may  only  be  found  in  experiments  with  high  levels  of 
background  noise,  where  local  disturbances  might  excite  its  development.  However, 
the  secondary  stability  analysis  of  Ragab  and  Wu  (1989)  shows  that  this  mode  can 
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become  important  when  compressibility  is  included,  and  may  be  important  down 
to  Mi  =  0.4. 


5.3  Effect  of  Mach  Number 

Most  of  the  three-dimensional  effects  observed  in  experiments  were  captured  in 
the  simulations  of  the  translative  mode  of  instability  in  the  previous  section.  It 
was  therefore  decided  to  use  this  combination  of  instability  waves  to  investigate  the 
effect  of  Mach  number.  Three  simulations  were  made,  at  Mach  numbers  Mi  =  0.4, 
0.8  and  1.05,  with  Reynolds  numbers  400,  600,  800  and  box  lengths  corresponding 
to  the  most  amplified  wavelengths  of  Lx  =  7.85,  13.37,  18.48  and  Lz  —  7.85,  13.37, 
12.465  respectively.  The  angles  of  the  oblique  instability  waves  were  45°  for  the 
M\  =  0.4  and  0.8  cases,  and  60°  for  the  case  Mi  =  1.05.  The  wave  combination  was 
the  fundamental  two-dimensional  wave  and  two  equal  and  opposite  oblique  waves 
with  the  same  wavelength.  The  amplitudes  were  chosen  to  be  Ai  =  Ai  —  0.025 
(equation  (5.1)),  with  <f>  =  n/2. 

The  effect  of  Mach  number  on  the  growth  of  the  mixing  layer  is  shown  on 
figure  5.18.  The  growth  in  the  energy  for  the  amplified  waves  is  shown  on  fig¬ 
ures  5.19  through  5.21.  At  Mi  =  0.4  it  can  be  seen  that  the  (1,0)  wave  is  the 
most  amplified  initially  and  is  always  more  energetic  than  the  oblique  waves.  At 
Mi  =  0.8  the  oblique  waves  (1,1)  and  (1,-1)  are  slightly  more  amplified  than 
the  two-dimensional  wave  and,  although  they  start  with  slightly  less  energy,  they 
soon  overtake  the  (1,0)  wave.  Again  the  linear  behavior  persists  in  the  non-linear 
regime  and  the  oblique  waves  have  a  higher  energy  content  in  the  developed  struc¬ 
ture  than  the  two-dimensional  waves,  although  at  this  intermediate  Mach  num¬ 
ber  both  two-dimensional  and  oblique  waves  are  important.  At  the  highest  Mach 
number,  Mi  =  1.05,  the  oblique  waves  have  a  much  larger  growth  rate  than  the 
two-dimensional  waves  and  by  the  end  of  the  simulation  the  energy  content  of  the 
(l,  1)  and  (1,  —  l)  modes  is  one  and  a  half  orders  of  magnitude  higher  than  the  (1,0) 
mode. 

In  each  of  the  simulations  the  linear  theory  correctly  predicts  the  initial  growth 
rate  of  the  instability  waves.  The  most  amplified  wave  is  also  the  most  important 
wave  in  the  developed  structure.  The  dominance  of  oblique  waves  at  high  Mach 
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numbers,  predicted  by  the  linear  stability  theory,  extends  to  the  non-linear  regime 
in  the  simulations  studied  here. 

To  illustrate  the  resolution  of  these  simulations  the  energy  content  of  all  the 
modes  is  shown  in  carpet  plots  of  E(kx,kz)  at  the  y  =  0  plane  (the  mixing  layer 
centerline)  on  figures  5.22  through  5.24  for  Mach  numbers  0.4,  0.8  and  1.05.  The 
drop  off  in  energy  from  the  most  energetic  wave  to  the  highest  wavenumbers  was 
kept  to  about  8  orders  of  magnitude  by  increasing  the  resolution  during  the  simu¬ 
lations.  When  the  highest  wavenumbers  became  important  the  simulations  had  to 
be  stopped.  The  plots  are  shown  for  the  last  time  step  in  the  simulation,  when  the 
resolution  was  considered  to  be  adequate. 

The  high  frequency  part  of  these  spectra  should  not  be  mistaken  for  small  eddies 
in  the  flow.  These  high  wavenumbers  arise  from  the  Fourier  representation  of  the 
large  structures  in  the  flow;  a  steep  gradient  at  some  point  in  the  structure  requires 
many  Fourier  modes  to  resolve  it.  For  example,  a  purely  two-dimensional  flow  would 
have  energy  only  along  the  kx  =  0  line.  At  low  Mach  number  the  carpet  plot  (figure 
5.22)  shows  oblique  ridges,  due  to  oblique  steep  gradients  that  need  to  be  resolved. 
The  plots  at  higher  Mach  numbers  (figures  5.23  and  5.24)  show  the  development 
of  ridges  in  the  kz  direction,  along  the  kx  =  0  line.  These  indicate  the  presence  of 
steep  spanwise  gradients  in  the  flow. 

Low  pressure  regions  are  associated  with  strong  rotation.  Perspective  views  of  a 
pressure  surface,  that  encloses  a  minima  of  pressure,  are  shown  on  figure  5.25  for 
M\  =  0.4  and  on  figures  5.27  and  5.28  for  Mach  numbers  0.8  and  1.05,  and  show  the 
change  in  large-scale  structure  in  the  flow  as  Mach  number  is  increased.  At  Mi  =  0.4 
we  have  the  translative  mode  discussed  in  the  previous  section,  with  the  vortex  tube 
oscillating  in  the  spanwise  direction.  Again  we  have  pre-collapse  streamwise  and 
vertical  vorticity  in  the  saddle  region  between  adjacent  rollers,  shown  by  the  surface 
of  streamwise  vorticity  on  figure  5.26. 

At  M\  =  0.8  (figure  5.27)  we  find  a  weakening  of  the  spanwise  structure,  and 
the  development  of  oblique  vortices  in  the  saddle  point  region  where  at  lower  Mach 
numbers  the  streamwise  braid  vortices  would  form.  By  Mi  =  1.05  (figure  5.28)  we 
find  that  the  two-dimensional  mode  has  all  but  vanished,  and  we  are  left  with  a 
pattern  of  four  regions  of  rotating  fluid.  There  is  one  pair  of  equal  and  opposite 
oblique  vortices  in  the  region  where  the  spanwise  vortex  was  located  at  lower  Mach 
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number  (x  =  0),  and  another  pair  of  counter-rotating  vortices,  with  opposite  sense 
to  the  first  two,  in  the  region  of  the  low  Mach  number  saddle  point  ( x  =  Lx/ 2). 

The  essential  change  in  structure  is  shown  on  figures  5.29  through  5.31,  where  two 
vortex  lLies  (lines  tangent  to  the  local  vorticity  vector)  per  structure  are  plotted.  At 
low  Mach  number  one  of  these  lines  is  the  spanwise  structure,  while  the  other  weaves 
back  and  forth  and  appears  as  the  streamwise  vortices  in  the  saddle  point  region. 
At  higher  Mach  number  both  vortex  lines  have  a  zig-zag  structure,  or  a  ‘double 
hairpin  with  peak-valley  splitting’  structure,  using  boundary-layer  terminology.  The 
pressure  minima  in  figures  5.27  and  5.28  show  the  positions  along  the  vortex  lines 
where  there  is  strong  rotation  taking  place.  At  other  regions,  for  example  near  the 
heads  of  the  hairpins,  there  is  very  little  rotation  of  fluid  taking  place.  The  presence 
of  a  vortex  line  in  this  region  should  not  be  mistaken  for  the  presence  of  a  vortex. 

The  sense  of  rotation  of  the  vortices  is  always  clockwise  if  a  cut  through  the  x  —  y 
plane  is  considered,  which  can  be  clearly  seen  in  the  mixture  fraction  field.  Cuts 
through  the  pressure  and  mixture  fraction  fields  at  the  plane  y  =  0  are  shown  on 
figure  5.32  for  M\  =  0.4  and  on  figures  5.33  and  5.34  for  Mi  =  0.8  and  M\  =  1.05 
respectively.  Fluid  from  below  the  mixing  layer  being  moved  upward  appears  as 
a  local  minima  in  the  mixture  fraction,  and  vice  versa  for  fluid  from  above  being 
moved  down.  To  enhance  the  effect  contour  plots  are  made  with  0.5  subtracted 
from  the  mixture  fraction,  which  lies  between  0  and  1.  Negative  values  are  shown 
with  dashed  contours,  labeling  fluid  from  below  the  mixing  layer  which  has  been 
moved  upwards.  Positive  values  are  shown  with  solid  contours  and  show  fluid  from 
above  that  has  been  moved  downwards. 

The  effect  of  Mach  number  on  the  bulging  mode  is  now  considered.  At  low  Mach 
numbers  this  mode  was  not  important,  while  at  high  Mach  numbers  the  structure 
is  expected  to  develop  from  oblique  instability  waves  alone,  with  no  effect  from  the 
relative  phase  of  any  two-dimensional  mode.  At  the  intermediate  Mach  number 
Mi  =  0.8  a  simulation  was  run  with  phase  <j>  =  0  (equation  (5.1)).  The  growth  in 
vorticity  thickness  and  mode  energies  for  this  simulation  are  shown  on  figures  5.35 
and  5.36,  and  the  developed  structure  on  figure  5.37.  As  with  the  other  simulation 
at  this  Mach  number  the  oblique  waves  are  the  most  amplified  instability  waves 
through  both  the  linear  and  non-linear  stages  of  roll-up.  The  final  structure,  though 
different  in  detail,  is  just  as  three-dimensional  as  that  which  developed  with  phase 
4>  =  7r/4.  It  also  contains  oblique  inclined  vortical  regions.  It  is  concluded  that  at 
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intermediate  Mach  numbers  the  phase  of  the  two-dimensional  wave  relative  to  the 
oblique  waves  has  an  effect  on  the  final  structure.  However  the  trend  towards  more 
three-dimensionality  at  higher  Mach  numbers  is  not  affected. 


5.4  Simulations  with  Random  Initial  Conditions 

The  purpose  of  these  simulations  is  to  check  that  the  linear  stability  theory 
is  correctly  predicting  the  most  amplified  waves  in  the  flow,  and  to  check  for  the 
presence  of  modes  other  than  those  found  in  the  linear  analysis.  The  box  length  was 
chosen  to  be  40  initial  vorticity  thicknesses  in  both  the  x  and  z  directions,  which  is 
long  enough  to  allow  approximately  5  of  the  most  amplified  waves  to  grow  at  Mach 
number  M\  =  0.4,  or  3  waves  at  =  0.8,  or  about  2  waves  at  the  highest  Mach 
number  M\  =  1.05.  The  initial  condition  was  specified  by  adding  a  small  random 
number  to  each  of  the  computational  variables  Q  —  (p,pu,pv,pw,e)  at  each  mesh 
point.  For  example 

p(x,y,z)  =  p(x,y,z)  +  Are  y  (5.3) 

where  r  is  a  random  number  uniformly  distributed  between  —0.5  and  0.5,  and  the 
amplitude  A  was  set  to  0.0001.  The  exponential  term  is  used  to  ensure  that  the 
disturbances  decay  to  zero  in  the  free-stream.  The  random  numbers  were  initially 
applied  to  a  16  x  33  x  16  grid,  which  was  then  extended  to  32  x  33  x  32  so  that 
random  numbers  in  the  highest  frequency  modes  could  not  alias  back  and  affect  the 
growth  rates  of  low  wavenumber  modes. 

The  simulations  were  run  for  the  same  conditions  as  in  the  previous  section 
i.e.  at  Mi  =  0.4  (Re  =  400),  Afx  =  0.8  (Re  =  600)  and  Mx  =  1.05  (Re  =  800). 
The  simulations  were  run  through  the  linear  regime  of  small  disturbance  growth. 
During  this  time  the  mean  flow  changes  by  viscous  diffusion.  This  is  illustrated 
by  the  variation  in  vorticity  thickness  during  these  simulations,  shown  on  figure 
5.38.  The  mean  flow  sets  the  length  scale  for  the  growth  of  small  disturbances, 
so  the  variation  in  the  mean  flow  means  that  different  wavelength  disturbances  are 
more  rapidly  amplified  at  different  times  during  the  simulations.  Specifically,  longer 
wavelengths  will  become  more  rapidly  amplified  as  the  simulations  proceed,  which 
needs  to  be  kept  in  mind  when  analyzing  the  results. 
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The  growth  in  energy  E(kx,  kz)  for  selected  modes  is  shown  on  figures  5.39  though 
5.41  for  the  Mach  numbers  M\  =  0.4,  0.8  and  1.05  respectively.  At  Mi  =  0.4  the 
growth  of  the  two-dimensional  modes  is  plotted.  Initially  the  (5,0)  mode  is  the  most 
amplified  of  those  shown,  in  agreement  with  linear  stability  theory.  However,  by 
the  end  of  the  simulation,  when  the  layer  has  grown  by  viscous  diffusion,  the  (3,0) 
mode  is  the  most  amplified.  Modes  with  wavelengths  longer  than  the  initially  most 
amplified  wavelength  have  growth  rates  that  increase  with  time,  whereas  modes  with 
shorter  wavelengths  have  growth  rates  that  diminish  with  time.  At  M\  =  0.8  the 
(3,4)  mode  is  more  amplified  than  the  (3,0),  (3,2)  for  (3,6)  modes.  At  M\  =  1.05 
the  (2,4)  mode  is  the  most  amplified  wave,  growing  more  strongly  than  the  (2,0), 
(2,2)  or  (2,6)  modes.  Both  these  results  are  in  general  agreement  with  the  linear 
stability  result  that  the  most  amplified  waves  satisfy  Mc  cos  0  =  0.6. 

The  final  flowfield  at  Mi  =  0.4  is  shown  on  figure  5.42  by  plotting  a  single  contour 
of  zero  vertical  velocity  in  the  x  —  z  plane  at  y  =  0.  This  contour  divides  fluid 
moving  upwards  from  fluid  moving  downwards  and  gives  an  idea  of  the  dominant 
structures  in  the  flow.  Clearly  there  is  a  preference  for  structures  oriented  in  the 
spanwise  direction,  though  there  is  no  strong  coherence. 

The  lack  of  the  strong  coherence  exhibited  here  at  Mi  =  0.4,  compared  to  incom¬ 
pressible  experiments,  may  have  many  causes.  In  the  simulations  all  waves  (0,0) 
through  (7,  7)  and  (7,  —7)  were  seeded,  so  there  are  many  waves  with  approximately 
the  correct  orientation  to  grow  nearly  as  strongly  as  the  most  amplified  wave.  The 
actual  amplitude  at  the  end  of  the  simulations  depends  upon  how  well  the  particular 
instability  mode  was  initialized.  Other  possible  causes  are  compressibility  effects, 
non-linear  effects,  or  the  variation  in  the  mean  flow  during  the  simulation.  It  may 
also  be  due  to  effects  peculiar  to  the  experiments.  For  example,  the  flow  is  espe¬ 
cially  receptive  to  disturbances  at  the  splitter  plate  edge,  which  is  two-dimensional. 
Also,  the  experiments  which  show  strong  spanwise  coherence  are  ‘clean’  experiments 
with  laminar  boundary-layers  coming  off  the  splitter  plate.  These  boundary-layers 
will  contain  instability  waves  which,  though  not  of  high  amplitude,  may  have  been 
growing  for  a  long  enough  time  to  have  sorted  out  the  two-dimensional  Tollmien- 
Schlichting  waves.  These  would  then  be  a  forcing  on  the  mixing  layer,  tending  to 
make  it  more  two-dimensional. 

The  changing  nature  of  the  most  unstable  waves  due  to  compressibility  is  clearly 
illustrated  in  the  sequence  of  figures  5.42  through  5.44.  These  show  the  v  =  0 
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contour  in  a  cut  through  the  x  —  z  plane  at  y  =  0  for  Mach  numbers  0.4,  0.8  and 
1.05  respectively.  At  M\  =  0.8  there  is  no  longer  any  tendency  towards  a  spanwise 
coherence,  and  waves  at  about  45°  are  most  common.  The  situation  is  clearer  at 
Afi  =  1.05  where  the  linear  theory  predicts  a  fairly  narrow  band  of  waves  around 
60°  to  be  most  amplified.  The  contours  of  v  =  0  on  figure  5.44  show  strong  evidence 
for  waves  at  this  angle.  There  are  regions  where  the  +  and  —  oblique  waves  are 
amplified  separately,  and  regions  where  both  are  amplified  at  the  same  region  in 
space. 

The  simulations  with  random  initial  conditions  confirmed  the  linear  stability 
finding  that  oblique  waves  become  the  most  amplified  waves  as  the  Mach  number 
of  the  mixing  layer  is  increased.  No  evidence  was  found  for  the  existence  of  modes 
other  than  those  already  considered  in  the  linear  stability  analysis. 

5.5  Structure  at  High  Mach  Number 

Using  results  from  the  previous  sections  we  can  make  some  predictions  about 
the  kind  of  large-scale  structures  which  may  grow  from  the  inflectional  instability 
in  the  mixing  layer  at  high  Mach  number,  especially  above  Mc  —  1  where  the 
oblique  waves  are  much  more  strongly  amplified  than  the  two-dimensional  waves. 
From  figure  5.44,  which  showed  the  structure  growing  from  low  amplitude  random 
noise,  we  see  that  there  are  regions  in  the  flow  where  one  of  the  oblique  waves 
seems  to  dominate  over  the  other,  and  regions  where  both  oblique  waves  seem  to 
be  about  equally  important.  Therefore  two  simulations  were  run,  one  with  a  single 
oblique  wave  and  the  other  with  a  pair  of  equal  and  opposite  oblique  waves.  The 
Mach  number  for  these  simulations  was  chosen  to  be  0.8,  so  that  a  lower  Reynolds 
number  of  400  could  be  used,  to  simulate  further  into  the  non-linear  development. 
The  two-dimensional  wave  was  not  included  in  the  initial  field,  so  these  structures 
are  expected  to  be  more  typical  of  flows  at  higher  Mach  numbers,  except  that  those 
structures  are  expected  to  become  more  oblique.  A  wave  angle  of  45°  was  selected, 
since  this  is  the  most  amplified  wave  at  Mi  =  0.8.  The  growth  in  vorticity  thickness 
is  shown  on  figure  5.45. 

The  developed  structure  arising  from  a  single  oblique  wave  is,  not  surprisingly, 
an  oblique  vortex.  The  structure  developing  from  two  equal  and  opposite  oblique 


waves  is  more  complex.  Pressure  surfaces  are  shown  on  figures  5.46  and  5.47  for  the 
two  cases.  In  the  structure  developing  from  the  pair  of  oblique  waves  there  are  four 
main  regions  of  rotating  fluid.  At  x  =  0  there  are  two  counter  rotating  vortices, 
inclined  in  y  relative  to  the  x  axis,  and  oblique  in  z  relative  to  the  x  axis.  There  are 
two  more  vortices  at  x  =  Lxj 2  similar  to  the  first  two,  but  with  the  opposite  sense 
of  u>x  rotation.  The  region  where  two  of  the  vortices  come  close  may  be  considered 
similar  to  the  hairpin  structures  found  in  transitional  boundary-layer  flows.  The 
induced  motion  of  the  head  of  the  hairpin  due  to  the  legs  is  alternately  up  or  down, 
explaining  the  inclined  nature  of  the  vortices.  In  fact,  the  vortices  become  more 
inclined  as  the  simulation  proceeds.  The  actual  rotation  at  the  heads  of  the  hairpins 
is  weak,  and  there  is  no  suggestion  of  rotation  of  the  head  of  one  hairpin  around 
the  tail  of  the  hairpin  beneath  it.  The  structure  can  be  thought  of  as  two  vortex 
lines,  passing  through  the  peaks  in  vorticity.  One  line  passes  through  the  vortices 
at  x  =  0  ,  z  =  Lz/A,3Lz/A  and  the  other  passes  through  the  vortices  at  x  =  Lxj 2, 
z  =  Z,z/4, 3Lz/4.  Perspective  and  top  views  of  the  vortex  lines  are  shown  on  figure 
5.48.  These  vortex  lines  are  staggered  in  the  streamwise  direction,  similar  to  the 
peak  valley  splitting  of  boundary-layer  transition.  However,  the  boundary  layer 
case  is  a  subharmonic  secondary  instability,  whereas  the  case  described  here  is  a 
fundamental  primary  instability. 

The  mixture  fraction  field  for  the  structure  resulting  from  two  equal  and  opposite 
oblique  waves  is  especially  rich  in  detail.  The  mixture  fraction,  /,  was  initially 
specified  by  a  hyperbolic  tangent  profile  (4.3),  and  tags  fluid  from  the  free-streams 
with  a  value  between  0  (lower  stream)  and  1  (upper  stream).  Contours  are  plotted 
of  /  —  0.5,  with  negative  contours  dashed,  so  that  solid  contours  mark  fluid  that 
originated  on  the  upper  side,  and  dashed  contours  mark  fluid  that  originated  on 
the  lower  side.  A  cut  through  the  x  —  z  plane  at  y  =  0  is  shown  on  figure  5.49.  The 
four  main  vortices  show  clearly.  The  regions  of  strong  mixture  fraction  gradient 
at  x  =  Lz/4,3Lz/4,  z  =  Lz/A^Lz/A  are  complex  three-dimensional  saddle  points, 
where  fluid  is  being  brought  to  rest,  and  high  pressure  ensues.  Cuts  of  mixture 
fraction  through  the  y  —  z  plane  at  x  —  Lx/A  and  x  =  Lxj 2  are  shown  on  figure 
5.50.  Mushroom  structures  are  found  at  x  =  0,  Lz/2  due  to  the  counter-rotating 
vortices.  Examples  of  cuts  in  the  x  —  y  plane  are  shown  on  figure  5.51  for  z  =  Lz/A 
and  z  =  Lz/2.  This  illustrates  that  care  is  needed  when  interpreting  experimental 
cuts  in  the  x  —  y  plane  or  span  averaged  photographs  such  as  Schlieren.  The  cuts 
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at  z  =  L2/4,3L2/4  are  at  first  glance  reminiscent  of  the  two-dimensional  low  Mach 
number  structure,  even  though  the  complete  structure  is  highly  three-dimensional. 
We  note  that  there  are  two  of  the  hairpin  structures  per  period,  which  may  explain 
the  reduction  in  structure  spacing  to  thickness  ratio  observed  by  Papamoschou 
[1986]  above  a  convective  Mach  number  of  0.8. 


5.6  Sensitivity  to  Initial  Conditions 

In  this  section  we  present  simulations  in  which  random  noise  was  added  to  the  ini¬ 
tial  condition.  The  objective  of  these  simulations  was  to  see  whether  the  structures 
presented  in  the  previous  sections  axe  strongly  affected  by  background  noise,  and  in 
particular  whether  non-linear  interactions  act  to  disrupt  the  organized  structure. 

A  simulation  at  M\  =  0.8  showed  that  the  growth  of  the  linear  eigenfunctions  was 
not  significantly  affected  by  the  presence  of  random  noise.  This  simulation  was  run 
with  the  same  Reynolds  number,  box  size  and  configuration  of  instability  waves  as 
in  section  5.3  (Re  —  600,  Lx  =  Lz  =  13.37,  forcing  with  a  two-dimensional  wave  and 
a  pair  of  equal  and  opposite  oblique  waves).  The  amplitude  of  each  wave  was  0.025. 
Random  noise  was  added  to  this  initial  condition  in  the  same  way  as  in  section  5.4, 
with  amplitude  A  =  0.025  (equation  (5.3)).  All  modes  from  (0,0)  through  (7,7)  and 
(7,-7)  were  seeded.  The  growth  in  mode  energies  of  the  unstable  waves  are  shown 
on  figure  5.52,  and  are  nearly  indistinguishable  from  the  growth  rates  shown  on 
figure  5.20,  where  there  was  no  random  noise.  The  final  structure  was  very  similar 
to  that  of  the  simulation  described  in  section  5.4  and  is  not  presented  again  here. 

A  more  stringent  test  was  run  in  which  the  initial  eigenfunctions  were  not  put 
into  the  simulations,  and  the  flow  had  to  sort  out  the  most  amplified  waves  from 
an  initial  field  consisting  only  of  random  noise.  The  same  Mach  number,  Reynolds 
number  and  box  size  were  used  as  in  section  5.3  and  in  the  preceding  paragraph.  The 
random  noise  was  added  with  amplitude  A  =  0.025,  seeding  modes  (0,0)  through 
(7,7)  and  (7,-7)  in  the  usual  way.  A  carpet  plot  of  the  initial  mode  energies  is  shown 
on  figure  5.53a,  and  a  cut  through  the  pressure  field  at  the  y  =  0  plane  is  shown  on 
figure  5.53 6.  The  simulation  was  run  forward  in  time  and  a  carpet  plot  and  pressure 
field  at  time  t  =  29.6  are  shown  on  figures  5.54a  and  5.546  respectively.  From  the 
carpet  plot  it  can  be  seen  that  the  energy  that  was  initially  put  into  modes  with 
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high  kx  has  decayed  very  rapidly.  The  peaks  that  have  appeared  in  the  carpet  plot 
are  the  unstable  waves.  From  linear  stability  theory  waves  (1,0),  (1,1),  (1,-1),  (1,2) 
and  (1,-2)  are  all  unstable.  Other  parts  of  the  carpet  plot  are  associated  either  with 
the  Fourier  representation  of  the  non-linear  structure  that  is  developing  from  the 
linear  instabilities,  or  with  the  decaying  initial  random  number  field. 

The  amplitude  of  each  wave  is  dependent  not  only  on  the  growth  rate,  but  also 
on  the  initial  random  field.  From  the  carpet  plot  (figure  5.54a)  it  appears  that 
the  (1,1)  mode  was  most  strongly  seeded,  followed  by  the  (1,-2)  and  (1,-1)  modes. 
As  time  advances  in  the  simulation  we  see  more  of  the  effect  of  amplification  rate, 
and  less  of  the  effect  of  the  initial  conditions.  The  pressure  field  shown  on  figure 
5.54ft  shows  that  one  of  the  oblique  waves  is  stronger  than  the  others  for  this  initial 
random  field. 

The  simulation  with  purely  random  initial  conditions  was  run  long  enough  for  the 
instability  waves  to  grow  to  near-saturation  conditions.  Organized  structure  was 
found.  At  time  t  =  52.0  the  carpet  plot  of  mode  energies,  and  the  pressure  field  at 
y  =  0  are  shown  on  figures  5.55a  and  5.55ft  respectively.  The  energy  spectra  has  now 
filled  out.  This  is  attributed  to  the  Fourier  modes  needed  to  resolve  the  large-scale 
structure  in  the  flowfield.  An  oblique  view  of  a  surface  of  constant  pressure  is  shown 
on  figure  5.56.  This  can  be  compared  with  figure  5.27  (structure  developing  from 
a  two-dimensional  wave  and  a  pair  of  oblique  waves  with  relative  phase  <f>  =  tt/2), 
figure  5.37  (same  but  with  (f>  =  0),  and  figure  5.47  (structure  developing  from  a  pair 
of  oblique  waves  alone).  The  developed  structure  is  similar  to  the  case  of  figure 
5.47,  where  the  structure  was  allowed  to  develop  from  a  pair  of  equal  and  opposite 
oblique  instability  waves. 

Figure  5.57  shows  how  the  linearly  unstable  waves  grow  from  out  of  the  random 
initial  condition.  There  is  an  adjustment  period,  extending  up  to  about  t  =  10, 
during  which  the  unstable  waves  emerge  from  the  background  noise.  Depending  on 
the  initial  field  different  waves  may  emerge  first.  It  is  apparent  that  the  (1,1)  mode 
was  most  strongly  seeded  in  this  case.  This  wave  grows  at  the  same  rate  as  the 
(1,-1)  wave  but  appears  to  start  from  an  earlier  time  (approximately  5  time  units). 
The  two-dimensional  wave  was  not  strongly  seeded  in  this  simulation,  and  does  not 
play  a  significant  role  in  the  evolution  of  this  flow. 

A  plot  of  the  growth  in  vorticity  thickness  (equation  (2.20))  is  shown  on  figure 
5.58a,  and  a  plot  of  the  growth  in  vorticity  thickness  based  on  a  mass-weighted 
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velocity  profile  is  shown  on  figure  5.58 b.  As  was  noted  in  section  4.3  the  vorticity 
thickness  is  very  sensitive  to  features  of  the  large-scale  structure.  Figure  5.58a 
shows  a  dip  in  the  vorticity  thickness  for  times  between  40  and  50,  while  figure 
5.58 b  shows  a  strong  increase  at  these  times.  These  effects  are  probably  specific  to 
this  simulation,  and  may  be  associated  with  the  ‘collapse’  of  vorticity  into  vortices, 
which  occurs  at  about  these  times. 

The  conclusions  of  these  simulations  are  (1)  that  the  linear  instability  processes 
are  not  very  sensitive  to  the  presence  of  background  noise,  and  (2)  that  an  organized 
structure  develops  from  simulations  started  with  random  initial  conditions.  This 
organized  structure  is  similar  to  that  which  was  presented  in  earlier  sections,  where 
the  structures  developed  from  combinations  of  the  instability  waves  with  the  highest 
linear  amplification  rates. 

The  apparent  ‘cleanliness’  of  these  simulations  is  not  due  to  the  choice  of  initial 
conditions,  but  to  the  dominance  of  the  primary  (inviscid  inflectional)  instability 
in  this  prototypical  free  shear  layer.  Even  with  random  initial  conditions  the  un¬ 
stable  modes  eventually  dominate  the  flow.  For  reasons  of  computational  expense 
the  simulations  have  not  been  run  to  later  times,  when  secondary  instabilities  will 
presumably  develop,  or  through  to  the  ‘mixing  transition’  which  (by  analogy  with 
the  incompressible  flow)  is  when  small  eddies  develop  in  the  mixing  layer.  However, 
even  beyond  these  stages  the  strong  instability  mechanism  will  certainly  persist,  and 
it  is  felt  that  the  large-scale  organized  structure  in  the  compressible  mixing  layer 
will  result  from  the  saturation  of  the  various  combinations  of  instability  waves,  as 
simulated  in  this  thesis. 


5.7  Chapter  Summary 

The  three-dimensional  simulations  presented  in  this  chapter  confirm  the  earlier 
linear  stability  finding  that  oblique  waves  are  more  amplified  than  two-dimensional 
waves  at  high  Mach  numbers.  Simulations  with  random  initial  conditions  showed 
that  the  structures  developing  from  linear  instabilities  become  more  oblique  as  Mach 
number  is  increased.  No  evidence  was  found  for  any  modes  of  instability  other  than 
those  already  found  in  the  linear  stability  analysis. 

At  convective  Mach  number  Mc  =  0.8  it  was  found  that  the  oblique  wave  was 
more  amplified  than  the  two-dimensional  wave  during  the  entire  simulation,  and 
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oblique  modes  contained  more  energy  than  two-dimensional  modes.  With  the  con¬ 
vective  Mach  number  above  1  (Mc  =  1.05)  it  was  found  that  the  two-dimensional 
instability  played  only  a  small  role  in  the  development  of  large-scale  structure. 

The  expected  structures  at  high  Mach  number  were  examined  in  more  detail  in 
two  more  simulations  -  one  of  a  single  oblique  wave  and  the  other  of  a  pair  of  equal 
and  opposite  oblique  waves.  The  oblique  wave  developed  into  an  oblique  vortex, 
while  the  pair  of  waves  developed  into  a  complex  structure  consisting  of  four  regions 
of  strong  rotation,  placed  along  two  vortex  lines  with  hairpin  shapes  in  the  spanwise 
direction,  and  with  peak-valley  splitting  in  the  streamwise  direction. 

Simulations  with  random  initial  conditions  were  run  to  check  the  sensitivity  of 
the  organized  structure  to  initial  conditions.  A  very  similar  structure  was  found  to 
develop.  The  growth  rate  of  the  linear  instabilities  was  found  to  be  insensitive  to 
the  presence  of  finite  amplitude  random  noise. 

No  shock  waves  were  found  in  any  of  the  three-dimensional  simulations,  even 
with  the  free-stream  Mach  number  above  1.  Whilst  this  does  not  prove  that  there 
will  never  be  shock  waves  in  mixing  layers  at  high  Mach  number,  it  does  appear 
that  the  flow  can  adjust  in  a  three-dimensional  manner  so  that  shock  waves  are  not 
required  at  high  Mach  numbers. 

The  expected  structure  of  the  compressible  mixing  layer  can  be  summarized  as 
follows.  Below  Mc  =  0.4  we  expect  the  usual  two-dimensional  rollers  that  have 
been  found  in  incompressible  experiments  and  simulations.  Above  Mc  =  1  we 
expect  to  find  the  double  hairpin  structure,  which  develops  from  pairs  of  oblique 
waves.  In  between  the  situation  is  more  complex,  since  a  broad  range  of  instability 
waves  are  strongly  amplified.  Between  Mc  =  0.4  and  Me  =  0.6  we  expect  that  the 
roller  structure  will  disappear.  Above  Mc  =  0.6  we  expect  to  see  more  and  more 
evidence  of  structures  with  strong  streamwise  vorticity,  such  as  might  develop  from 
combinations  of  two-dimensional  and  oblique  instability  waves. 
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CHAPTER  0 


Conclusions  and  Recommendations 


This  work  has  been  concerned  with  a  numerical  study  of  the  compressible  mix¬ 
ing  layer.  The  methods  used  were  linear  stability  analysis  and  direct  numerical 
simulation  of  the  compressible  Navier-Stokes  equations.  The  stability  equations 
were  solved  by  a  shooting  method.  The  full  equations  were  solved  by  an  explicit 
code,  with  derivatives  evaluated  spectrally  or  with  high  order  finite  differences. 
The  conclusions  are  divided  into  the  three  main  areas  of  study:  linear  stability, 
two-dimensional  simulations  and  three-dimensional  simulations. 


Linear  Stability  Theory: 

•  Linear  stability  theory  can  be  used  to  predict  mixing  layer  growth  rates.  The 
relation  61  ~  Icqlmax  gives  the  correct  trends  in  growth  rate  due  to  velocity 
ratio,  density  ratio  and  Mach  number.  The  amplification  rate  of  the  most  am¬ 
plified  wave  has  to  be  computed  using  spatial  stability  theory,  with  a  solution 
to  the  boundary-layer  equations  as  the  base  state. 

•  The  experimental  finding  of  Papamoschou  [1989]  that  convective  velocities  do 
not  match  the  prediction  of  the  Uc  formula  (1.2)  was  reproduced  in  the  linear 
theory,  using  the  phase  speed  of  the  neutral  modes  as  a  prediction  of  convective 
velocities  of  large-scale  structures. 

•  Oblique  waves  were  found  to  be  more  amplified  above  a  convective  Mach 
number  of  0.6.  A  simple  relation,  Afccos0  =  0.6  was  proposed  to  give  an 
approximate  orientation  of  the  most  amplified  waves  in  the  flow. 

•  Several  different  methods  of  achieving  high  convective  Mach  numbers  were 
compared.  There  was  no  overall  collapse  of  growth  rates  at  high  Mach  num¬ 
bers,  indicating  that  the  convective  Mach  number  Mc  =  {U[  —  £/|)/(cJ  +  c2) 
may  be  only  a  first  order  measure  of  compressibility  effects. 
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Two-Dimensional  Simulations: 


•  The  growth  rate  of  the  two-dimensional  mixing  layer  was  observed  to  drop 
rapidly  as  the  Mach  number  was  increased,  nonlinearly  as  well  as  linearly. 

•  The  reduction  in  growth  rate  was  associated  with  a  change  in  shape  of  the 
developed  vortices.  These  became  more  elongated  in  the  streamwise  direction 
as  Mach  number  was  increased. 

•  Shock  waves  were  observed  in  two-dimensional  simulations  above  a  convective 
Mach  number  of  0.7. 

•  Peak  strain  rate  in  the  simulations  was  found  to  be  always  the  same  order  as 
the  global  strain  rate,  indicating  that  the  global  strain  rate  can  be  used  to 
predict  the  magnitude  of  local  strain  rate  in  the  mixing  layer. 

•  Simulations  with  a  density  ratio  of  0.2  were  performed  for  a  case  where  the 
Uc  formula  (1.2)  would  predict  zero  convection  speed.  It  was  found  that  the 
structures  did  indeed  move,  in  agreement  with  the  phase  speed  of  the  neutrally 
stable  mode. 

•  The  supersonic  modes  of  instability,  which  are  the  only  unstable  modes  in  two- 
dimensions  at  high  Mach  number,  were  simulated.  The  modes  were  confirmed 
to  be  radiating,  with  a  pattern  of  shock  and  expansion  waves  forming  on  the 
side  of  the  mixing  layer  relative  to  which  the  instability  was  supersonic.  These 
modes  were  found  to  be  vortical,  and  did  lead  to  roll-up.  However,  the  growth 
rate  was  very  small. 


Three-Dimensional  Simulations: 

•  At  low  Mach  numbers  it  was  found  that  the  phase  of  a  two-dimensional  wave 
relative  to  a  pair  of  equal  and  opposite  oblique  waves  could  be  chosen  to  give 
the  ‘bulging’  or  ‘translative’  secondary  stability  modes  of  Pierrehumbert  and 
Widnall  [1982].  Simulations  of  these  cases  confirmed  their  finding  that  there 
was  significant  instability  only  for  the  translative  mode. 

•  Raising  the  Mach  number,  with  the  phasing  chosen  to  give  the  translative 
mode,  was  found  to  lead  to  a  change  in  structure.  The  linear  stability  result 
that  oblique  waves  are  more  amplified  than  two-dimensional  waves  was  found 
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to  carry  over  into  the  nonlinear  regime.  Above  a  convective  Mach  number  of 
0.6  the  oblique  modes  of  the  developed  structure  were  found  to  contain  more 
energy  than  the  two-dimensional  mode. 

•  Above  a  convective  Mach  number  of  1  there  was  found  to  be  very  little  influ¬ 
ence  of  the  two-dimensional  wave  on  the  developed  structure. 

•  No  shock  waves  were  found  in  the  three-dimensional  simulations,  even  above 
a  convective  Mach  number  of  1. 

•  Simulations  with  random  initial  conditions  confirmed  the  linear  stability  find¬ 
ing  of  oblique  waves  being  more  amplified  at  high  Mach  numbers.  The  linear 
theory  was  found  to  predict  well  the  angle  of  the  most  unstable  modes  as 
Mach  number  was  increased.  No  evidence  was  found  for  any  modes  of  insta¬ 
bility  other  than  those  already  found  in  the  linear  stability  analysis. 

•  The  addition  of  random  noise  did  not  significantly  change  the  growth  rate  of 
linear  instability  waves. 

•  Typical  structures  which  may  be  found  in  the  mixing  layer  at  high  Mach 
numbers  were  computed,  based  on  the  nonlinear  evolution  of  oblique  instability 
waves.  A  single  oblique  wave  led  to  an  oblique  vortex.  A  pair  of  equal  and 
opposite  oblique  waves  led  to  a  structure  with  four  regions  of  strong  rotation, 
arranged  along  two  kinked  spanwise  vortex  tubes,  staggered  in  the  streamwise 
direction.  The  structure  resembles  a  peak-valley  splitting  arrangement  of  the 
hairpin  structures  found  in  boundary-layer  transition. 

•  A  simulation  starting  with  random  initial  conditions  developed  a  large-scale 
structure  very  similar  to  that  computed  from  the  non-linear  development  of 
combinations  of  instability  waves. 

•  The  mixing  layer  below  a  convective  Mach  number  Mc  =  0.4  is  expected  to 
have  a  ‘roller’  structure  similar  to  the  incompressible  flow.  Between  Mc  =  0.4 
and  Mc  =  1.0  we  expect  to  find  a  complex  structure,  due  to  many  instability 
waves  being  equally  amplified,  but  with  a  tendency  to  stronger  streamwise 
and  oblique  vortices  as  the  Mach  number  is  raised.  Above  Mc  =  1  we  expect 
to  find  a  cleaner  structure  of  pairs  of  hairpins,  which  develop  from  pairs  of 
oblique  instability  waves. 
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Open  Questions  for  Future  Research 

Numerical  simulation  of  compressible  turbulent  flows  is  a  relatively  new  area  of 
research  and  many  questions  remain  to  be  answered.  Future  work  could  be  directed 
towards  some  of  the  following  questions,  which  arose  from  this  work. 

Linear  Instability  Model: 

•  What  are  the  limitations  of  the  model  discussed  in  section  2.3? 

The  immediate  need  is  for  more  experimental  measurements  of  growth  rates 
and  convective  velocities  as  a  function  of  density  ratio  and  Mach  number.  Sta¬ 
bility  analysis  of  actual  self-similar  experimental  velocity  and  density  profiles 
may  help  to  isolate  some  of  the  limitations  of  the  simple  model. 

Subharmonic  Growth  Mechanism: 

•  What  is  the  high  Mach  number  analog  of  the  vortex  pairing  process? 

There  is  presumed  to  be  some  mechanism  by  which  the  largest  structures  in 
flow  merge  so  that  the  mixing  layer  cross-stream  length  scale  can  grow  and  the 
mixing  layer  can  be  self-similar.  At  high  Mach  numbers  this  probably  involves 
oblique  subharmonic  waves.  It  is  not  clear  what  kind  of  ‘oblique  pairing’  of 
vortices,  or  ‘merging  of  hairpin  structures’,  will  take  place.  Simulation  of  this 
process  will  be  expensive,  since  twice  the  resolution  in  each  direction  will  be 
required,  compared  to  the  current  work. 

Mixing  Transition: 

•  How  do  small  eddies  form  in  a  free-shear  layer? 

The  mechanism  is  not  understood  even  in  incompressible  flow,  but  the  experi¬ 
mental  observation  is  of  a  large  increase  in  formation  of  products  of  a  chemical 
reaction  at  some  point  in  the  mixing  layer  development  (Breidenthal  [1981]). 
The  resolution  required  for  the  Reynolds  numbers  at  which  mixing  transition 
occurs  is  probably  near  current  supercomputer  capability. 

Statistics: 

•  How  do  we  make  engineering  computations  of  compressible  turbulent  flows? 
Traditional  models  require  a  knowledge  of  the  time-averaged  character  of  the 
flow.  Statistical  information  from  compressible  turbulence  simulations  could 
be  used  to  point  to  the  relevant  terms  in  the  equations  which  need  to  be 
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modelled  differently  at  higher  Mach  numbers.  However  it  is  felt  that  no 
simulation  has  yet  produced  a  statistically  self-similar  mixing  layer,  even  in 
two-dimensional  incompressible  flow.  A  large  enough  sample  of  the  largest 
structures  ( e.g .  60-100)  would  have  to  be  accumulated  to  get  reliable  statistics 
(Sandham  and  Reynolds  [1989]).  A  logical  sequence  for  progress  on  statistics 
seems  to  be:  two-dimensional  incompressible,  three-dimensional  incompress¬ 
ible,  and  finally  three-dimensional  compressible.  Compressible  simulations 
typically  require  twice  the  storage  and  three  times  the  CPU  time  of  incom¬ 
pressible  simulations. 

Mixing: 

•  Does  the  scalar  pdf  change  at  the  higher  Mach  numbers? 

The  change  in  large-scale  structure  suggests  that  the  scalar  probability  density 
function  (pdf)  will  be  different  at  higher  Mach  numbers.  A  change  in  the  pdf 
also  means  that  fast  chemical  reactions  will  behave  differently  to  the  low  Mach 
number  flow.  Scalar  pdf’s  can  be  accumulated  from  the  simulations  presented 
here.  However,  they  are  of  limited  value  since  we  are  not  simulating  through 
the  mixing  transition,  and  have  no  small-scales  doing  the  mixing.  Also  we 
do  not  have  a  self-similar  mixing  layer,  so  the  pdf  varies  with  time  in  the 
simulations. 

Chemical  Reactions: 

•  How  to  make  supersonic  combustion  more  efficient? 

We  need  to  understand  how  heat  release  affects  the  compressible  mixing  layer 
flow.  One  possible  avenue  for  research  is  to  follow  the  procedure  of  this  thesis. 
First  of  all,  the  stability  problem  could  be  solved  for  the  compressible  mixing 
layer,  with  an  artificial  temperature  profile  to  reflect  the  effects  of  heat  release. 
Then  direct  numerical  simulations  could  be  performed  with  two  species,  in¬ 
cluding  reaction  rate  terms.  The  effect  of  changes  in  the  instability  of  the  flow 
could  be  studied,  as  well  as  the  effect  of  finite-rate  chemistry. 
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A  Direct  Method  for  Linear  Stability  Analysis 

In  Chapter  2  the  linearized  equations  were  reduced  to  a  single  equation  which 
was  then  solved  by  shooting.  In  this  appendix  an  alternative  procedure  is  de¬ 
scribed,  which  makes  use  of  a  spectral  representation  of  the  mean  and  perturbed 
flow  to  formulate  the  linear  stability  problem  as  a  matrix  eigenvalue  problem.  This 
method  was  developed  in  collaboration  with  J.  H.  Chen;  results  for  compressible 
wake  instability  can  be  found  in  Chen  tt  al.  [1989]. 

We  begin  with  the  linearized  equations  which  were  derived  in  section  2.1.2.  As 
before  the  notation  D  is  used  for  the  d/dy  operator.  The  linearized  continuity 
equation  is: 

pi{au  —  w)  +  vDp  +  p[t(au  +  0u>)  +  Dv]  =  0  (A.l) 

After  substituting  for  the  linearized  perfect  gas  law,  equation  (2.39),  the  linearized 
momentum  equations  in  the  x,  y  and  z  directions  become: 


_r ./  _  ^ ^  —ia(pT  +  pT) 

p[t(au  -  w)u  -f-  vDu\  - - ^ - 

(A.2) 

,,  -{pDT  +  TDp  +  pDT  +  TDp) 

")»  - 

(A.3) 

-  >  -  -i0(pT  +  pT) 

=  lM2 

(A-  4) 

and  the  linearized  energy  equation  is: 

p[t(au  -  u)T  -|-  vDT]  —  —{7  -  l)[*(au  +  0w)  +  Dv]  (A. 5) 

The  first  step  is  to  rearrange  these  equations  into  the  form  of  a  matrix  multiplied 
by  a  vector,  with  the  desired  eigenvalue  isolated.  For  temporal  stability  calculations 
we  wish  to  isolate  u/,  so  we  write  the  above  equations  as: 

A(f  =  wi  (-4-6) 
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where  x  is  a  column  vector  of  the  eigenfunctions  x  =  (p,u,t),u/,T)r  and  the  matrix 
At  is  given  by 
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For  spatial  stability  we  isolate  a,  and  write: 

Aax  =  aBax 

where  the  matrix  As  is  given  by: 


{A.  8) 
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and  the  matrix  is  given  by: 
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(A. 10) 


The  next  step  is  to  derive  a  matrix  form  of  the  operator  D.  The  mapped  Fourier 
method  of  Cain  et  al.  [1984]  was  selected  since  it  allows  solution  over  an  infinite 
physical  domain.  We  start  with  the  definition  of  a  Fourier  transform  pair: 

4  =  (4.11) 

y= o 
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N/  2—1 

<%)  =  E 

k=-N/2 


(A.  12) 


where  Af  is  the  number  of  points,  rjj  are  the  grid  points,  <j>  is  the  function  and  <j>  is 
the  transformed  function.  Cain  et  al.  [1984]  defined  a  cotangent  mapping: 


y  =  — accotan(27r^) 


(A. 13) 


where  ac  is  a  constant  that  determines  the  amount  of  stretching  in  the  mapping. 
Under  this  mapping  the  mixing  layer  physical  domain  —  oo  <  y  <  oo  maps  onto 
0  <  rj  <  1/2.  The  computational  domain  0  <  rj  <  1  is  periodic,  allowing  the  use  of 
Fourier  methods. 

Derivatives  are  calculated  as: 


d(f>  1  d(f> 
dy  h'  dr) 

where  h!  =  dy/drj  is  the  metric.  For  Cain’s  mapping  this  is  given  by: 


(A. 14) 


1  1 


h!  27ra, 


— *  sin2(27r?7)  =  — —  fl 

27rac  47rac  L 


ei2(2nrj)  +  e~i2(2irrj) 


(A.15) 


Substituting  this  into  (A. 14),  using  (A. 12)  and  combining  terms,  gives  an  expression 
for  the  derivative  at  grid  point  j: 


dy 


N/2—1 


i{k-2)~  i{k  +  2)2 


H  [ikh  -  _  o 

1  2ac  ,  .  L  2  2 


k=—N/2+l 


h+2  ]eik2^  (A.  16) 


where  all  terms  with  wavenumbers  outside  the  summation  limits  are  discarded. 

The  next  step  is  to  substitute  in  for  0  from  equation  (A. 11),  and  compare  with 
a  matrix  formulation: 

t  =  dlj<f>{Vj)  U-17) 
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The  d\j  are  terms  of  a  matrix,  given  by: 


ffl 2—1 

dtj  =  — Y  ft'/c  -  t(fe~2)et4^  -  »(fc  +  2)e— «4^y]e»2irM»»<-»Ty)  (>U8) 
2acTV  ^  L  2  2  J  v  ' 

ib=-J\r/2+l 

This  expression  can  be  simplified  using  trigonometry  and  geometric  series. 

For  l  =  j  we  obtain: 

d[j  =  -  ^  N-  sin(4?r ijy)  (A.19a) 


and  for  l  ^  j: 

dlj  _ 1  ,-.c°s<4/’';;).(-1)--/cot(^i))  _  {A.m) 

h  4  ac  K  1  \  N  J  acN  sin[^Vi] 

As  a  last  step  we  use  the  symmetry  of  the  Cain  mapping  to  reduce  the  size  of 
the  matrix  from  N  x  N  to  N*  x  N*,  where  N*  =  N/2  +  1.  A  new  matrix  e/y  is 
defined  as  follows. 

For  j  =  1  and  j  =  JV*: 

elj  =  (A. 20a) 

and  for  j  =  2  to  N*  —  1: 

elj  =  dlj  +  dl(N+2-j)  (A. 20b) 

This  matrix  is  used  to  represent  derivative  terms  in  matrices  At  (equation  (A.7)) 
and  A8  (equation  (A.9)).  The  result  is  a  5 N*  x  5 N*  matrix  (5x5  matrices  at 
each  grid  point).  The  eigenvalues  and  eigenvectors  of  this  matrix  were  found  using 
standard  (IMSL)  subroutines. 

The  performance  of  the  method  was  checked  by  running  the  test  case  of  table  2.2. 
Curves  of  amplification  rate  against  wavenumber  are  shown  on  figure  A.l  for  several 
values  of  N*t  and  compared  with  the  result  from  the  shooting  method  of  Chapter  2. 
It  can  be  seen  that  convergence  of  the  direct  method  is  good  for  the  most  amplified 
wave,  but  very  poor  near  the  neutral  mode  and  at  very  low  u>.  In  general,  it  was 
found  that  weakly  amplified  waves  were  difficult  to  resolve  with  the  direct  method. 
At  high  Mach  numbers  it  was  difficult  to  capture  the  weakly  amplified  instability 
modes,  even  using  TV*  as  high  as  60. 
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Figure  1.2  Time-developing  mixing  layer. 
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Figure  2.2  Velocity  and  density  profiles  computed  from  the  boundary-layer  equa¬ 
tions,  showing  the  effect  of  density  ratio:  (a)  velocity  profiles,  (b) 
density  profiles. 


Figure  2.3  Spatial  amplification  rate  plotted  against  A  =  (£/*  —  U^JIU  r  +  ffa). 

showing  the  combined  effects  of  velocity  and  density  ratio:  (a)  profiles 
from  boundary-layer  equations,  (b)  hyperbolic  tangent. 
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Figure  2.4  Variation  in  maximum  spatial  amplification  rate  with  density  ratio 
compared  with  models  for  growth  rate  from  Brown  [1974]  and  Dimo- 
takis  [1986],  for  velocity  ratios:  (a)  U2  =  0,  (b)  U2  =  0.5. 


(a) 


(b) 


Figure  2.5  Plot  of  amplification  rate  against  the  functional  form  of  the  models 
of  (a)  Brown  [1974]  (model  1)  and  (b)  Dimotakis  [1986]  (model  2). 
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Figure  2.10  Variation  of  temporal  amplification  rate  with  angle  of  disturbance  for 
the  hyperbolic  tangent  velocity  profile. 


Figure  2.11  Variation  of  spatial  amplification  rate  with  angle  of  disturbance  for 
the  mixing  layer  with  Tjj  =  1.0  and  U2  =  0.5,  with  velocity  profile 
from  solution  of  the  boundary-layer  equations. 
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Figure  2.12  Variation  of  the  amplification  rate  of  the  most  unstable  wave  with 
Mach  number,  (a)  comparing  the  most  amplified  wave  with  the  most 
amplified  two-dimensional  wave,  (b)  extension  to  Mc  =  3.2  for  the 
most  unstable  wave. 
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Figure  2.13  Variation  of  spatial  amplification  rate  of  the  most  unstable  wave  with 
convective  Mach  number  for  three  different  methods  of  varying  con¬ 
vective  Mach  number. 
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Figure  2.14  Growth  rate  data  from  Papamoschou  [1986]  normalized  by  the  am¬ 
plification  rate  of  the  most  unstable  wave,  for  the  given  velocity  and 
density  ratio,  at  zero  Mach  number 
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Figure  2.15  Variation  in  amplification  rate  with  Mach  number  M\  for  the  subsonic 
mode,  the  ‘fast’  supersonic  mode,  and  the  ‘slow’  supersonic  mode. 
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Figure  2.19  Eigenfunctions  of  pressure  at  Mi  =  2.2  for  (a)  the  fast  supersonic 
mode,  (b)  the  slow  supersonic  modes  and  (c)  the  subsonic  mode. 


(e)  max=1.20e-3  min=-1.20e-3 


(f)  max=1.37e-3  min=-1.41e-3 


(g)  max=4.68e-4  min=-4.78e-4 


Figure  2.22  Contour  plots  from  linear  eigenfunctions  at  M\  =  0.6  (a)  loz,  (b) 
vz/Pi  (c)  density,  (d)  pressure,  (e)  dilatation,  (f)  dilatation  term,  (g) 
baroclinic  term. 


Figure  2.23  Schematic  of  the  s 
developing  mixing 


Figure  2.24  Amplification  rate  versus  frequency,  defining  the  most  amplified  wave, 
and  the  neutral  wave. 
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Figure  2.25  Plot  of  Mc\  versus  Mc 2  showing  a  comparison  between  the  theoretical 
Ue  from  equation  (1.2),  the  linear  stability  prediction  (the  phase  speed 
of  the  neutral  mode),  and  the  experimental  data  from  Papamoschou 
[1989]. 
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Decay  of  u  velocity  at  x  =  0.5  z  =  0.25,  for  Taylor-Green  problem. 


(a)  max=0.500  min=-0.500 


(b)  max=1.59  min=-6.95e-3 


(a)  max=0.500  min=-0.500 


(b)  max— 1.59  min— -0.226 


Figure  3.6  Developed  structure  for  Ly  —  6  (a)  mixture  fraction  (b)  vorticity. 
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Figure  3.7  Energy  spectrum  for  the  case  Ly  =  10,  M\  =  0.4. 
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Figure  4.2  Effect  of  Reynolds  number  on  the  growth  history  of  vorticity  thickness 
at  Mi  =  0.4. 


Figure  4.6  Long-time  behavior  of  vorticity  thickness  at  M\  =  0.4. 
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Figure  4.7  Structure  at  long  time  (a)  mixture  fraction,  (b)  vorticity. 
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Figure  4.9 
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Comparison  of  vorticity  thickness  measures,  based  on  mean  velocity 
profile  (R),  or  mass-weighted  mean  velocity  profile  (MW). 


Effect  of  Mach  n  umber  on  the  growth  in  mode  energy  E  of  the  most 
amplified  disturbance. 
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(a)  max=0.500  min =-0.500 


(b)  max=18.3  min=16.8 


(c)  max=1.81  min=-4.44e-3 


(d)  max =1.91  min=-4.44e-3 


Figure  4.10  Developed  structure  at  Mi 
(c)  uz,  (d) 

=  0.2  (a)  mixture  fraction,  (b)  pressure, 
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(c)  max=1.59  min=-6.95e-3 


(d)  max =1.93  min=-6.96e-3 


Figure  4.11  Developed  structure  at  M\  =  0.4  (a)  mixture  fraction,  (b)  pressure, 
(c)  w„  (d)  w,/p. 
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Figure  4.12  Developed  structure  at  M\ 
(c)  w„  (d)  uig/p. 


0.6  (a)  mixture  fraction,  (b)  pressure, 
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(a)  max=0.500  min=-0.500 


(b)  max =1.66  min=0.698 


X  X 


(c)  max=1.08  min=-2.82e-2  (d)  max=1.74  min=-2.54e-2 


Figure  4.13  Developed  structure  at  M\  =  0.8  (a)  mixture  fraction,  (b)  pressure, 
(c)  loz,  (d)  uz/p. 
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(b)  max =0.107  min=-9.36e-2 


Figure  4.14  Developed  structure  at  Mi  =  0.6  (a)  density,  (b)  dilatation,  (c)  di- 
latational  term  in  vorticity  equation,  (d)  baroclinic  term. 
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Figure  4.15  Choice  of  phase  between  subharmonic  and  fundamental  modes. 
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Figure  4.17  Growth  in  mode  energy  E,  comparing  M\  =  0.2  with  M\  =  0.6. 
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(c)  max=1.76  min=-2.89e-4  (d)  max=1.82  min=-2.89e-4 
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Figure  4.18  Step  1  in  pairing  process  at  Mj  =  0.2  (a)  mixture  fraction,  (b)  pres¬ 
sure,  (c)  (d)  ujx/p. 
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Figure  4.19  Step  2  in  pairing  process  at  M i  =  0.2  (a)  mixture  fraction,  (b)  pres¬ 
sure,  (c)  u/z,  (d) 
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(a)  max=0.500  min— -0.500 


(b)  max=18.2  min=16.2 


(c)  max=1.54  min=-4.89e-2  (d)  max=1.65  min=-4.88e-2 


Figure  4.20  Step  3  in  pairing  process  at  M\  =  0.2  (a)  mixture  fraction,  (b)  pres¬ 
sure,  (c)  u»„  (d)  ug/p. 


(c)  max=1.48  min=-7.81e-2  (d)  max=1.58  min=-7.83e-2 


Figure  4.21  Step  4  in  pairing  process  at  M\  =  0.2  (a)  mixture  fraction,  (b)  pres¬ 
sure,  (c)  (d)  w,/p. 
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Figure  4.22  Step  5  in  pairing  process  at  M\  —  0.2  (a)  mixture  fraction,  (b)  pres¬ 
sure,  (c)  u2,  (d)  u)zjp. 
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Figure  4.26  Step  4  in  pairing  process  at  =  0.6  (a)  mixture  fraction 


Figure  4.27  Comparison  of  vorticity  thickness  growth  at  two  Mach  numbers  where 
the  initial  wavelengths  are  the  same. 
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Figure  4.29  Plot  of  strain  rate  and  mixture  fraction  at  Mi  =  0.6,  t  =  18.2  (a 
mixture  fraction,  (b)  strain  rate  S. 


Figure  4.30  Plot  of  strain  rate  and  mixture  fraction  at  Mi  =  0.6,  t  =  24.0  (a 
mixture  fraction,  (b)  strain  rate  5. 
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Figure  4.34  Growth  in  mode  energy  E,  for  temperature  ratios  1  and  2. 
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Q  Generation  of  clockwise  vorticity 


Generation  of  counter-clockwise  vorticity 


Figure  4.36  Generation  of  baroclinic  torques  in  the  saddle  region. 
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I  Figure  4.37  Developed  structure  at  Mi  =  0.6,  T2  =  5  (a)  mixture  fraction,  (b) 

pressure,  (c)  uZi  (d)  w*/p. 

I 

I 
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Figure  4.41  Growth  of  vorticity  thickness  for  supersonic  instability  mode. 


I 

Figure  4.42  Growth  of  mode  energy  for  supersonic  instability  mode. 
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Figure  4.43  Non-linear  structure  developing  from  supersonic  mode  instability:  (a) 
mixture  fraction,  (b)  vorticity,  (c)  density,  (d)  pressure. 


Figure  4.44  Pressure  contours  showing  development  of  shock  waves  from  super¬ 
sonic  mode  instability. 
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Figure  5.1 
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Figure  5.2 


Linear  temporal  amplification  rate  versus  wave  angle  at  Mach  num¬ 
bers  0.95  and  1.05. 


Effect  of  Reynolds  number  on  the  temporal  amplification  rate  of 
oblique  disturbances  in  the  compressible  mixing  layer. 
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(a)  max =4. 51  min=4.42 


(a)  max =4. 50  min=4.42 


(b)  max=0.258  min=-0.258 


Figure  5.4  Initial  condition  for  the  translative  mode  of  instability  (a)  pressure 
minimum  at  x  =  0,  maximum  at  Lz/2  (b)  ux  (c)  uy  (d)  uz. 


Figure  5.5  Vorticity  thickness  growth  for  the  bulging  and  translative  modes  at 
Mi  =  0.4. 


Figure  5.6  Growth  in  energy  for  the  bulging  instability  at  Mi  =  0.4.  (1,0)  is  the 
two-dimensional  wave.  (1,1)  and  (1,-1)  are  the  oblique  waves. 


Figure  5.7  Growth  in  energy  for  the  translative  instability  at  Mi  =  0.4.  (1,0)  is 
the  two-dimensional  wave.  (1,1)  and  (1,-1)  are  the  oblique  waves. 


181 


Figure  5.13  Vorticity  thickness  growth  for  the  helical  subharmonic  mode  at  M\ 
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Figure  5.14  Growth  in  energy  for  the  bulging  instability  at  A/j  =  0.4.  (2,0)  is  the 
two-dimensional  wave.  (1,1)  and  (1,-1)  are  the  oblique  waves. 
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of  the  pressure  minima  at  time  t  =  12.63  for  the 
lie  mode  of  instability,  showing  vortex  cores  at  sat- 
idamental  instability. 
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Figure  5.16  Perspective  view  of  the  pressure  minima  at  time  t  =  27.90  for  the  he¬ 
lical  subharmonic  mode  of  instability,  showing  the  final  subharmonic 
structure. 
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(a)  max=0.489  min=-0.489 


(b)  max— 4.98  min=3.59 


Figure  5.17  Cuts  at  y  =  0  through  the  final  structure  developed  from  the  helical 
subharmonic  mode  (a)  mixture  fraction  (b)  pressure. 
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Figure  5.19  Growth  in  mode  energies  at  =  0.4  (1,0)  is  the  two-dimensional 
mode.  (1,1)  and  (1,-1)  are  the  oblique  waves. 
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Figure  5.20  Growth  in  mode  energies  at  Mi  =  0.8  (1,0)  is  the  two-dimensional 
mode.  (1,1)  and  (1,-1)  are  the  oblique  waves. 
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Figure  5.21  Growth  in  mode  energies  at  M\  =  1.05  (1,0)  is  the  two-dimensional 
mode.  (1,1)  and  (1,-1)  are  the  oblique  waves. 


Figure  5.25  Surface  of  constant  pressure  showing  the  rotational  region  in  the  struc¬ 
ture  developing  at  Afj  =  0.4. 


Figure  5.26  Perspective  view  of  streamwise  vorticity  in  the  structure  that  devel¬ 
oped  at  Mi  =  0.4 
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Figure  5.27  Surface  of  constant  pressure  showing  the  rotational  region  in  the  struc¬ 


ture  developing  at  M\  =  0.8. 
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the  peaks  of  vorticity  for  the  developed  structure 
>ective  view  (b)  top  view. 
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(a)  max=0.464  min=-0.464 


(b)  max=1.43  min =0.651 


Figure  5.33  Cuts  at  y  =  0  through  the  final  structure  developed  at  My  =  0.8  (a) 
mixture  fraction  (b)  pressure. 
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Figure  5.34  Cuts  at  y  =  0  through  the  final  structure  developed  at  M\  =  1.05  (a) 
mixture  fraction  (b)  pressure. 
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Figure  5.35  Vorticity  thickness  growth  for  simulation  of  the  bulging  mode  at  M\  = 
0.8. 


Figure  5.36  Growth  in  energy  for  the  bulging  instability  mode  at  M\  =  0.8. 
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Figure  5.37  Surface  of  constant  pressure  showing  the  rotational  region  in  the  struc 
ture  developing  from  the  bulging  mode  at  M\  =  0.8. 


Figure  5.39  Growth  in  energy  for  selected  modes  from  the  simulation  beginning 
with  random  numbers  at  Mj  =  0.4 


Figure  5.42  Flowfield  at  the  end  of  the  linear  stage  of  instability  growth  at  M\ 
0.4,  shown  by  the  v  =  0  contour. 


Figure  5.43  Flowfield  at  the  end  of  the  linear  stage  of  instability  growth  at  Mi 
0.8,  shown  by  the  v  =  0  contour. 
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Figure  5.44  Flowfield  at  the  end  of  the  linear  stage  of  instability  growth  at  Mi 
1.05,  shown  by  the  v  =  0  contour. 
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Figure  5.45  Vorticity  thickness  growth  for  a  single  versus  a  pair  of  oblique  waves. 
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Figure  5.52  Growth  in  mode  energies  of  the  unstable  waves  in  the  present 
background  noise. 
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Figure  5.53  Simulation  with  random  initial  conditions  at  time  t  =  0  (a)  carpet 
plot  of  mode  energies,  (b)  pressure  cut  at  y  =  0. 
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Figure  5.54  Simulation  with  random  initial  conditions  at  time  t  =  29.' 
plot  of  mode  energies,  (b)  pressure  cut  at  y  =  0. 
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Figure  5.56  Surface  of  constant  pressure  for  the  simulation  beginning  with  purely 
random  initial  conditions. 
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Figure  5.57  Growth  in  mode  energies  of  the  unstable  waves  in  the  simulation  be¬ 
ginning  with  purely  random  initial  conditions. 
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Figure  A.l  Illustration  of  the  convergence  of  the  direct  method.  The  number  of 
points  used  ( N *)  is  given  in  the  legend.  The  correct  solution,  as  found 
by  the  shooting  method,  is  also  shown. 
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